13 #include <Eigen/Eigen>
38 namespace SmallDeformation
40 template <
typename ShapeFunction,
41 typename IntegrationMethod,
43 SmallDeformationLocalAssemblerMatrixNearFracture<ShapeFunction,
46 SmallDeformationLocalAssemblerMatrixNearFracture(
48 std::size_t
const n_variables,
50 std::vector<unsigned>
const& dofIndex_to_localIndex,
51 bool const is_axially_symmetric,
52 unsigned const integration_order,
55 n_variables * ShapeFunction::NPOINTS * DisplacementDim,
56 dofIndex_to_localIndex),
57 _process_data(process_data),
58 _integration_method(integration_order),
60 _is_axially_symmetric(is_axially_symmetric)
66 DisplacementDim>(e, is_axially_symmetric,
69 unsigned const n_integration_points =
72 _ip_data.reserve(n_integration_points);
78 for (
unsigned ip = 0; ip < n_integration_points; ip++)
80 _ip_data.emplace_back(solid_material);
82 auto const& sm = shape_matrices[ip];
84 ip_data.dNdx = sm.dNdx;
85 ip_data.integration_weight =
87 sm.integralMeasure * sm.detJ;
90 static const int kelvin_vector_size =
92 ip_data._sigma.setZero(kelvin_vector_size);
93 ip_data._eps.setZero(kelvin_vector_size);
96 ip_data._sigma_prev.resize(kelvin_vector_size);
97 ip_data._eps_prev.resize(kelvin_vector_size);
99 ip_data._C.resize(kelvin_vector_size, kelvin_vector_size);
116 template <
typename ShapeFunction,
typename IntegrationMethod,
119 ShapeFunction, IntegrationMethod,
120 DisplacementDim>::assembleWithJacobian(
double const t,
double const dt,
121 Eigen::VectorXd
const& local_u,
122 Eigen::VectorXd& local_b,
123 Eigen::MatrixXd& local_J)
125 assert(_element.getDimension() == DisplacementDim);
127 auto const N_DOF_PER_VAR = ShapeFunction::NPOINTS * DisplacementDim;
128 auto const n_fractures = _fracture_props.size();
129 auto const n_junctions = _junction_props.size();
130 auto const n_enrich_var = n_fractures + n_junctions;
132 using BlockVectorType =
133 typename Eigen::VectorXd::FixedSegmentReturnType<N_DOF_PER_VAR>::Type;
134 using BlockMatrixType =
135 Eigen::Block<Eigen::MatrixXd, N_DOF_PER_VAR, N_DOF_PER_VAR>;
152 auto local_b_u = local_b.segment<N_DOF_PER_VAR>(0);
153 std::vector<BlockVectorType> vec_local_b_g;
154 for (
unsigned i = 0; i < n_enrich_var; i++)
156 vec_local_b_g.push_back(
157 local_b.segment<N_DOF_PER_VAR>(N_DOF_PER_VAR * (i + 1)));
160 auto local_J_uu = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(0, 0);
161 std::vector<BlockMatrixType> vec_local_J_ug;
162 std::vector<BlockMatrixType> vec_local_J_gu;
163 std::vector<std::vector<BlockMatrixType>> vec_local_J_gg(n_enrich_var);
164 for (
unsigned i = 0; i < n_enrich_var; i++)
166 auto sub_ug = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
167 0, N_DOF_PER_VAR * (i + 1));
168 vec_local_J_ug.push_back(sub_ug);
170 auto sub_gu = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
171 N_DOF_PER_VAR * (i + 1), 0);
172 vec_local_J_gu.push_back(sub_gu);
174 for (
unsigned j = 0; j < n_enrich_var; j++)
176 auto sub_gg = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
177 N_DOF_PER_VAR * (i + 1), N_DOF_PER_VAR * (j + 1));
178 vec_local_J_gg[i].push_back(sub_gg);
182 auto const nodal_u = local_u.segment<N_DOF_PER_VAR>(0);
183 std::vector<BlockVectorType> vec_nodal_g;
184 for (
unsigned i = 0; i < n_enrich_var; i++)
186 auto sub =
const_cast<Eigen::VectorXd&
>(local_u).segment<N_DOF_PER_VAR>(
187 N_DOF_PER_VAR * (i + 1));
188 vec_nodal_g.push_back(sub);
194 unsigned const n_integration_points =
195 _integration_method.getNumberOfPoints();
202 for (
unsigned ip = 0; ip < n_integration_points; ip++)
206 auto& ip_data = _ip_data[ip];
207 auto const& w = _ip_data[ip].integration_weight;
209 auto const& N = ip_data.N;
210 auto const& dNdx = ip_data.dNdx;
213 Eigen::Vector3d
const ip_physical_coords(
215 std::vector<double>
const levelsets(
217 _fracID_to_local, ip_physical_coords));
222 for (
unsigned i = 0; i < n_enrich_var; i++)
224 nodal_total_u += levelsets[i] * vec_nodal_g[i];
228 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
232 ShapeFunction::NPOINTS,
234 dNdx, N, x_coord, _is_axially_symmetric);
237 auto const& eps_prev = ip_data._eps_prev;
238 auto const& sigma_prev = ip_data._sigma_prev;
240 auto& eps = ip_data._eps;
241 auto& sigma = ip_data._sigma;
242 auto& state = ip_data._material_state_variables;
244 eps.noalias() = B * nodal_total_u;
246 variables[
static_cast<int>(
254 variables_prev[
static_cast<int>(
258 variables_prev[
static_cast<int>(
260 .emplace<double>(_process_data._reference_temperature);
262 auto&& solution = _ip_data[ip]._solid_material.integrateStress(
263 variables_prev, variables, t, x_position, dt, *state);
267 OGS_FATAL(
"Computation of local constitutive relation failed.");
271 std::tie(sigma, state, C) = std::move(*solution);
275 local_b_u.noalias() -= B.transpose() * sigma * w;
276 for (
unsigned i = 0; i < n_enrich_var; i++)
278 vec_local_b_g[i].noalias() -=
279 levelsets[i] * B.transpose() * sigma * w;
283 local_J_uu.noalias() += B.transpose() * C * B * w;
285 for (
unsigned i = 0; i < n_enrich_var; i++)
288 vec_local_J_ug[i].noalias() +=
289 B.transpose() * C * (levelsets[i] * B) * w;
292 vec_local_J_gu[i].noalias() +=
293 (levelsets[i] * B.transpose()) * C * B * w;
295 for (
unsigned j = 0; j < n_enrich_var; j++)
298 vec_local_J_gg[i][j].noalias() +=
299 (levelsets[i] * B.transpose()) * C * (levelsets[j] * B) * w;
305 template <
typename ShapeFunction,
306 typename IntegrationMethod,
311 computeSecondaryVariableConcreteWithVector(
312 double const , Eigen::VectorXd
const& )
315 const int n = DisplacementDim == 2 ? 4 : 6;
316 Eigen::VectorXd ele_stress = Eigen::VectorXd::Zero(n);
317 Eigen::VectorXd ele_strain = Eigen::VectorXd::Zero(n);
319 unsigned const n_integration_points =
320 _integration_method.getNumberOfPoints();
321 for (
unsigned ip = 0; ip < n_integration_points; ip++)
323 auto& ip_data = _ip_data[ip];
325 ele_stress += ip_data._sigma;
326 ele_strain += ip_data._eps;
328 ele_stress /= n_integration_points;
329 ele_strain /= n_integration_points;
331 (*_process_data._mesh_prop_stress_xx)[_element.getID()] = ele_stress[0];
332 (*_process_data._mesh_prop_stress_yy)[_element.getID()] = ele_stress[1];
333 (*_process_data._mesh_prop_stress_zz)[_element.getID()] = ele_stress[2];
334 (*_process_data._mesh_prop_stress_xy)[_element.getID()] = ele_stress[3];
335 if (DisplacementDim == 3)
337 (*_process_data._mesh_prop_stress_yz)[_element.getID()] = ele_stress[4];
338 (*_process_data._mesh_prop_stress_xz)[_element.getID()] = ele_stress[5];
341 (*_process_data._mesh_prop_strain_xx)[_element.getID()] = ele_strain[0];
342 (*_process_data._mesh_prop_strain_yy)[_element.getID()] = ele_strain[1];
343 (*_process_data._mesh_prop_strain_zz)[_element.getID()] = ele_strain[2];
344 (*_process_data._mesh_prop_strain_xy)[_element.getID()] = ele_strain[3];
345 if (DisplacementDim == 3)
347 (*_process_data._mesh_prop_strain_yz)[_element.getID()] = ele_strain[4];
348 (*_process_data._mesh_prop_strain_xz)[_element.getID()] = ele_strain[5];
Definition of the Node class.
Definition of the Point3d class.
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)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
MechanicsBase< DisplacementDim > & selectSolidConstitutiveRelation(std::map< int, std::unique_ptr< MechanicsBase< DisplacementDim >>> const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
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)
MathLib::Point3d computePhysicalCoordinates(MeshLib::Element const &e, Eigen::MatrixBase< Derived > const &shape)
std::vector< double > uGlobalEnrichments(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)
BMatrixType computeBMatrix(DNDX_Type const &dNdx, N_Type const &N, const double radius, const bool is_axially_symmetric)
Fills a B-matrix based on given shape function dN/dx values.
Coordinates mapping matrices at particular location.
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N