13 #include <Eigen/Eigen>
25 namespace SmallDeformation
27 template <
typename ShapeFunction,
typename IntegrationMethod,
29 SmallDeformationLocalAssemblerFracture<ShapeFunction, IntegrationMethod,
31 SmallDeformationLocalAssemblerFracture(
33 std::size_t
const n_variables,
35 std::vector<unsigned>
const& dofIndex_to_localIndex,
36 bool const is_axially_symmetric,
37 unsigned const integration_order,
40 n_variables * ShapeFunction::NPOINTS * DisplacementDim,
41 dofIndex_to_localIndex),
42 _process_data(process_data),
43 _integration_method(integration_order),
46 DisplacementDim>(e, is_axially_symmetric,
47 _integration_method)),
52 unsigned const n_integration_points =
55 _ip_data.reserve(n_integration_points);
59 auto frac_id =
_process_data._map_materialID_to_fractureID[mat_id];
74 for (
unsigned ip = 0; ip < n_integration_points; ip++)
81 ip_data.integration_weight =
83 sm.integralMeasure * sm.detJ;
84 ip_data.h_matrices.setZero(DisplacementDim,
85 ShapeFunction::NPOINTS * DisplacementDim);
92 ip_data.w.setZero(DisplacementDim);
93 ip_data.sigma.setZero(DisplacementDim);
96 ip_data.sigma_prev.resize(DisplacementDim);
97 ip_data.w_prev.resize(DisplacementDim);
99 ip_data.C.resize(DisplacementDim, DisplacementDim);
102 ip_data.aperture_prev = ip_data.aperture0;
108 template <
typename ShapeFunction,
typename IntegrationMethod,
111 ShapeFunction, IntegrationMethod,
112 DisplacementDim>::assembleWithJacobian(
double const t,
double const ,
113 Eigen::VectorXd
const& local_u,
114 Eigen::VectorXd& local_b,
115 Eigen::MatrixXd& local_J)
117 auto const N_DOF_PER_VAR = ShapeFunction::NPOINTS * DisplacementDim;
118 auto const n_fractures = _fracture_props.size();
119 auto const n_junctions = _junction_props.size();
120 auto const n_enrich_var = n_fractures + n_junctions;
134 using BlockVectorType =
135 typename Eigen::VectorXd::FixedSegmentReturnType<N_DOF_PER_VAR>::Type;
136 using BlockMatrixType =
137 Eigen::Block<Eigen::MatrixXd, N_DOF_PER_VAR, N_DOF_PER_VAR>;
139 std::vector<BlockVectorType> vec_local_b_g;
140 for (
unsigned i = 0; i < n_enrich_var; i++)
142 vec_local_b_g.push_back(
143 local_b.segment<N_DOF_PER_VAR>(N_DOF_PER_VAR * i));
145 std::vector<std::vector<BlockMatrixType>> vec_local_J_gg(n_enrich_var);
146 for (
unsigned i = 0; i < n_enrich_var; i++)
148 for (
unsigned j = 0; j < n_enrich_var; j++)
150 auto sub_gg = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
151 N_DOF_PER_VAR * i, N_DOF_PER_VAR * j);
152 vec_local_J_gg[i].push_back(sub_gg);
156 std::vector<Eigen::VectorXd> vec_nodal_g;
157 for (
unsigned i = 0; i < n_enrich_var; i++)
159 auto sub =
const_cast<Eigen::VectorXd&
>(local_u).segment<N_DOF_PER_VAR>(
161 vec_nodal_g.push_back(sub);
169 int const index_normal = DisplacementDim - 1;
170 auto const& R = _fracture_property->R;
172 unsigned const n_integration_points =
173 _integration_method.getNumberOfPoints();
178 for (
unsigned ip = 0; ip < n_integration_points; ip++)
182 auto& ip_data = _ip_data[ip];
183 auto const& integration_weight = ip_data.integration_weight;
184 auto const& H = ip_data.h_matrices;
185 auto& mat = ip_data.fracture_material;
186 auto& sigma = ip_data.sigma;
187 auto const& sigma_prev = ip_data.sigma_prev;
189 auto const& w_prev = ip_data.w_prev;
191 auto& state = *ip_data.material_state_variables;
192 auto& N = _secondary_data.N[ip];
194 Eigen::Vector3d
const ip_physical_coords(
197 _fracture_property->fracture_id, _fracture_props, _junction_props,
198 _fracID_to_local, ip_physical_coords));
202 Eigen::VectorXd nodal_gap(N_DOF_PER_VAR);
204 for (
unsigned i = 0; i < n_enrich_var; i++)
206 nodal_gap += levelsets[i] * vec_nodal_g[i];
210 w.noalias() = R * H * nodal_gap;
213 ip_data.aperture = ip_data.aperture0 + w[index_normal];
216 mat.computeConstitutiveRelation(
217 t, x_position, ip_data.aperture0,
218 Eigen::Matrix<double, DisplacementDim, 1>::Zero(),
222 w_prev, w, sigma_prev, sigma, C, state);
225 for (
unsigned i = 0; i < n_enrich_var; i++)
227 vec_local_b_g[i].noalias() -= levelsets[i] * H.transpose() *
228 R.transpose() * sigma *
233 for (
unsigned i = 0; i < n_enrich_var; i++)
235 for (
unsigned j = 0; j < n_enrich_var; j++)
238 vec_local_J_gg[i][j].noalias() +=
239 (levelsets[i] * H.transpose() * R.transpose()) * C *
240 (levelsets[j] * R * H) * integration_weight;
246 template <
typename ShapeFunction,
typename IntegrationMethod,
250 computeSecondaryVariableConcreteWithVector(
const double t,
251 Eigen::VectorXd
const& local_u)
253 auto const N_DOF_PER_VAR = ShapeFunction::NPOINTS * DisplacementDim;
254 auto const n_fractures = _fracture_props.size();
255 auto const n_junctions = _junction_props.size();
256 auto const n_enrich_var = n_fractures + n_junctions;
258 auto const& R = _fracture_property->R;
262 int const index_normal = DisplacementDim - 1;
264 unsigned const n_integration_points =
265 _integration_method.getNumberOfPoints();
270 std::vector<Eigen::VectorXd> vec_nodal_g;
271 for (
unsigned i = 0; i < n_enrich_var; i++)
273 auto sub =
const_cast<Eigen::VectorXd&
>(local_u).segment<N_DOF_PER_VAR>(
275 vec_nodal_g.push_back(sub);
278 for (
unsigned ip = 0; ip < n_integration_points; ip++)
282 auto& ip_data = _ip_data[ip];
283 auto const& H = ip_data.h_matrices;
284 auto& mat = ip_data.fracture_material;
285 auto& sigma = ip_data.sigma;
286 auto const& sigma_prev = ip_data.sigma_prev;
288 auto const& w_prev = ip_data.w_prev;
290 auto& state = *ip_data.material_state_variables;
291 auto& b_m = ip_data.aperture;
292 auto& N = _secondary_data.N[ip];
294 Eigen::Vector3d
const ip_physical_coords(
297 _fracture_property->fracture_id, _fracture_props, _junction_props,
298 _fracID_to_local, ip_physical_coords));
302 Eigen::VectorXd nodal_gap(N_DOF_PER_VAR);
304 for (
unsigned i = 0; i < n_enrich_var; i++)
306 nodal_gap += levelsets[i] * vec_nodal_g[i];
310 w.noalias() = R * H * nodal_gap;
313 b_m = ip_data.aperture0 + w[index_normal];
317 "Element {:d}, gp {:d}: Fracture aperture is {:g}, but it must "
319 _element.getID(), ip, b_m);
323 mat.computeConstitutiveRelation(
324 t, x_position, ip_data.aperture0,
325 Eigen::Matrix<double, DisplacementDim, 1>::Zero(),
329 w_prev, w, sigma_prev, sigma, C, state);
334 HMatricesType::ForceVectorType::Zero(DisplacementDim);
336 HMatricesType::ForceVectorType::Zero(DisplacementDim);
338 for (
unsigned ip = 0; ip < n_integration_points; ip++)
340 auto& ip_data = _ip_data[ip];
341 ele_b += ip_data.aperture;
343 ele_sigma += ip_data.sigma;
345 ele_b /= n_integration_points;
346 ele_w /= n_integration_points;
347 ele_sigma /= n_integration_points;
348 (*_process_data._mesh_prop_b)[_element.getID()] = ele_b;
349 (*_process_data._mesh_prop_w_n)[_element.getID()] = ele_w[index_normal];
350 (*_process_data._mesh_prop_w_s)[_element.getID()] = ele_w[0];
351 (*_process_data._mesh_prop_fracture_stress_normal)[_element.getID()] =
352 ele_sigma[index_normal];
353 (*_process_data._mesh_prop_fracture_stress_shear)[_element.getID()] =
virtual constexpr unsigned getDimension() const =0
Get dimension of the mesh element.
virtual std::size_t getID() const final
Returns the ID of the element.
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
VectorType< DisplacementDim > ForceVectorType
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
std::vector< double > duGlobalEnrichments(std::size_t this_frac_id, std::vector< FractureProperty * > const &frac_props, std::vector< JunctionProperty * > const &junction_props, std::unordered_map< int, int > const &fracID_to_local, Eigen::Vector3d const &x)
MathLib::Point3d computePhysicalCoordinates(MeshLib::Element const &e, Eigen::MatrixBase< Derived > const &shape)
void computeHMatrix(N_Type const &N, HMatrixType &H)
Fills a H-matrix based on given shape function.
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N
ParameterLib::Parameter< double > const & aperture0
Initial aperture.