14#include <range/v3/range/conversion.hpp>
15#include <range/v3/view/transform.hpp>
39namespace SmallDeformation
41template <
typename ShapeFunction,
44SmallDeformationLocalAssemblerMatrixNearFracture<ShapeFunction,
47 SmallDeformationLocalAssemblerMatrixNearFracture(
49 std::size_t
const n_variables,
51 std::vector<unsigned>
const& dofIndex_to_localIndex,
53 bool const is_axially_symmetric,
56 n_variables * ShapeFunction::NPOINTS * DisplacementDim,
57 dofIndex_to_localIndex),
58 _process_data(process_data),
59 _integration_method(integration_method),
61 _is_axially_symmetric(is_axially_symmetric)
67 DisplacementDim>(e, is_axially_symmetric,
70 unsigned const n_integration_points =
73 _ip_data.reserve(n_integration_points);
79 for (
unsigned ip = 0; ip < n_integration_points; ip++)
81 _ip_data.emplace_back(solid_material);
83 auto const& sm = shape_matrices[ip];
85 ip_data.dNdx = sm.dNdx;
86 ip_data.integration_weight =
88 sm.integralMeasure * sm.detJ;
91 static const int kelvin_vector_size =
93 ip_data._sigma.setZero(kelvin_vector_size);
94 ip_data._eps.setZero(kelvin_vector_size);
97 ip_data._sigma_prev.resize(kelvin_vector_size);
98 ip_data._eps_prev.resize(kelvin_vector_size);
100 ip_data._C.resize(kelvin_vector_size, kelvin_vector_size);
112 ranges::views::transform(
115 ranges::to<std::vector>;
118template <
typename ShapeFunction,
int DisplacementDim>
121 DisplacementDim>::assembleWithJacobian(
double const t,
double const dt,
122 Eigen::VectorXd
const& local_u,
123 Eigen::VectorXd& local_b,
124 Eigen::MatrixXd& local_J)
126 assert(_element.getDimension() == DisplacementDim);
128 auto const N_DOF_PER_VAR = ShapeFunction::NPOINTS * DisplacementDim;
129 auto const n_fractures = _fracture_props.size();
130 auto const n_junctions = _junction_props.size();
131 auto const n_enrich_var = n_fractures + n_junctions;
133 using BlockVectorType =
134 typename Eigen::VectorXd::FixedSegmentReturnType<N_DOF_PER_VAR>::Type;
135 using BlockMatrixType =
136 Eigen::Block<Eigen::MatrixXd, N_DOF_PER_VAR, N_DOF_PER_VAR>;
153 auto local_b_u = local_b.segment<N_DOF_PER_VAR>(0);
154 std::vector<BlockVectorType> vec_local_b_g;
155 for (
unsigned i = 0; i < n_enrich_var; i++)
157 vec_local_b_g.push_back(
158 local_b.segment<N_DOF_PER_VAR>(N_DOF_PER_VAR * (i + 1)));
161 auto local_J_uu = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(0, 0);
162 std::vector<BlockMatrixType> vec_local_J_ug;
163 std::vector<BlockMatrixType> vec_local_J_gu;
164 std::vector<std::vector<BlockMatrixType>> vec_local_J_gg(n_enrich_var);
165 for (
unsigned i = 0; i < n_enrich_var; i++)
167 auto sub_ug = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
168 0, N_DOF_PER_VAR * (i + 1));
169 vec_local_J_ug.push_back(sub_ug);
171 auto sub_gu = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
172 N_DOF_PER_VAR * (i + 1), 0);
173 vec_local_J_gu.push_back(sub_gu);
175 for (
unsigned j = 0; j < n_enrich_var; j++)
177 auto sub_gg = local_J.block<N_DOF_PER_VAR, N_DOF_PER_VAR>(
178 N_DOF_PER_VAR * (i + 1), N_DOF_PER_VAR * (j + 1));
179 vec_local_J_gg[i].push_back(sub_gg);
183 auto const nodal_u = local_u.segment<N_DOF_PER_VAR>(0);
184 std::vector<BlockVectorType> vec_nodal_g;
185 for (
unsigned i = 0; i < n_enrich_var; i++)
187 auto sub =
const_cast<Eigen::VectorXd&
>(local_u).segment<N_DOF_PER_VAR>(
188 N_DOF_PER_VAR * (i + 1));
189 vec_nodal_g.push_back(sub);
195 unsigned const n_integration_points =
196 _integration_method.getNumberOfPoints();
203 for (
unsigned ip = 0; ip < n_integration_points; ip++)
207 auto& ip_data = _ip_data[ip];
208 auto const& w = _ip_data[ip].integration_weight;
210 auto const& N = ip_data.N;
211 auto const& dNdx = ip_data.dNdx;
214 Eigen::Vector3d
const ip_physical_coords(
216 std::vector<double>
const levelsets(
218 _fracID_to_local, ip_physical_coords));
223 for (
unsigned i = 0; i < n_enrich_var; i++)
225 nodal_total_u += levelsets[i] * vec_nodal_g[i];
229 NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
233 ShapeFunction::NPOINTS,
235 dNdx, N, x_coord, _is_axially_symmetric);
238 auto const& eps_prev = ip_data._eps_prev;
239 auto const& sigma_prev = ip_data._sigma_prev;
241 auto& eps = ip_data._eps;
242 auto& sigma = ip_data._sigma;
243 auto& state = ip_data._material_state_variables;
245 eps.noalias() = B * nodal_total_u;
257 variables_prev.
temperature = _process_data._reference_temperature;
259 auto&& solution = _ip_data[ip]._solid_material.integrateStress(
260 variables_prev, variables, t, x_position, dt, *state);
264 OGS_FATAL(
"Computation of local constitutive relation failed.");
268 std::tie(sigma, state, C) = std::move(*solution);
272 local_b_u.noalias() -= B.transpose() * sigma * w;
273 for (
unsigned i = 0; i < n_enrich_var; i++)
275 vec_local_b_g[i].noalias() -=
276 levelsets[i] * B.transpose() * sigma * w;
280 local_J_uu.noalias() += B.transpose() * C * B * w;
282 for (
unsigned i = 0; i < n_enrich_var; i++)
285 vec_local_J_ug[i].noalias() +=
286 B.transpose() * C * (levelsets[i] * B) * w;
289 vec_local_J_gu[i].noalias() +=
290 (levelsets[i] * B.transpose()) * C * B * w;
292 for (
unsigned j = 0; j < n_enrich_var; j++)
295 vec_local_J_gg[i][j].noalias() +=
296 (levelsets[i] * B.transpose()) * C * (levelsets[j] * B) * w;
302template <
typename ShapeFunction,
308 computeSecondaryVariableConcreteWithVector(
309 double const , Eigen::VectorXd
const& )
312 const int n = DisplacementDim == 2 ? 4 : 6;
313 Eigen::VectorXd ele_stress = Eigen::VectorXd::Zero(n);
314 Eigen::VectorXd ele_strain = Eigen::VectorXd::Zero(n);
316 unsigned const n_integration_points =
317 _integration_method.getNumberOfPoints();
318 for (
unsigned ip = 0; ip < n_integration_points; ip++)
320 auto const& ip_data = _ip_data[ip];
322 ele_stress += ip_data._sigma;
323 ele_strain += ip_data._eps;
325 ele_stress /= n_integration_points;
326 ele_strain /= n_integration_points;
328 (*_process_data._mesh_prop_stress_xx)[_element.getID()] = ele_stress[0];
329 (*_process_data._mesh_prop_stress_yy)[_element.getID()] = ele_stress[1];
330 (*_process_data._mesh_prop_stress_zz)[_element.getID()] = ele_stress[2];
331 (*_process_data._mesh_prop_stress_xy)[_element.getID()] = ele_stress[3];
332 if (DisplacementDim == 3)
334 (*_process_data._mesh_prop_stress_yz)[_element.getID()] = ele_stress[4];
335 (*_process_data._mesh_prop_stress_xz)[_element.getID()] = ele_stress[5];
338 (*_process_data._mesh_prop_strain_xx)[_element.getID()] = ele_strain[0];
339 (*_process_data._mesh_prop_strain_yy)[_element.getID()] = ele_strain[1];
340 (*_process_data._mesh_prop_strain_zz)[_element.getID()] = ele_strain[2];
341 (*_process_data._mesh_prop_strain_xy)[_element.getID()] = ele_strain[3];
342 if (DisplacementDim == 3)
344 (*_process_data._mesh_prop_strain_yz)[_element.getID()] = ele_strain[4];
345 (*_process_data._mesh_prop_strain_xz)[_element.getID()] = ele_strain[5];
Definition of the Node class.
Definition of the Point3d class.
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > mechanical_strain
std::variant< std::monostate, Eigen::Matrix< double, 4, 1 >, Eigen::Matrix< double, 6, 1 > > stress
virtual std::size_t getID() const final
Returns the ID of the element.
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
unsigned getNumberOfPoints() const
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
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)
Eigen::Vector3d 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