39namespace SmallDeformationNonlocal
45template <
typename ShapeMatrixType>
48 std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>>
N;
51template <
typename ShapeFunction,
int DisplacementDim>
81 bool const is_axially_symmetric,
88 unsigned const n_integration_points =
91 _ip_data.reserve(n_integration_points);
94 auto const shape_matrices =
96 DisplacementDim>(e, is_axially_symmetric,
99 auto& solid_material =
104 auto* ehlers_solid_material =
dynamic_cast<
107 if (ehlers_solid_material ==
nullptr)
110 "The SmallDeformationNonlocalLocal process supports only "
111 "Ehlers material at the moment. For other materials the "
112 "interface must be extended first.");
115 for (
unsigned ip = 0; ip < n_integration_points; ip++)
117 _ip_data.emplace_back(*ehlers_solid_material);
119 auto const& sm = shape_matrices[ip];
122 sm.integralMeasure * sm.detJ;
125 ip_data.dNdx = sm.dNdx;
128 ip_data.sigma.setZero(
135 ip_data.sigma_prev.resize(
138 ip_data.eps_prev.resize(
149 double const* values,
150 int const integration_order)
override
152 if (integration_order !=
156 "Setting integration point initial conditions; The integration "
157 "order of the local assembler for element {:d} is different "
158 "from the integration order in the initial condition.",
167 if (name ==
"kappa_d")
177 std::string
const& name, std::vector<double>
const& value)
override
179 if (name ==
"kappa_d_ip")
181 if (value.size() != 1)
184 "CellData for kappa_d initial conditions has wrong number "
185 "of components. 1 expected, got {:d}.",
194 double const internal_length2 =
_process_data.internal_length_squared;
195 return (distance2 > internal_length2)
197 : (1 - distance2 / (internal_length2)) *
198 (1 - distance2 / (internal_length2));
205 DisplacementDim>>>
const& local_assemblers)
override
210 unsigned const n_integration_points =
213 std::vector<double> distances;
219 for (
unsigned k = 0; k < n_integration_points; k++)
225 auto const& xyz =
_ip_data[k].coordinates;
228 for (
auto const search_element_id : search_element_ids)
230 auto const& la = local_assemblers[search_element_id];
231 la->getIntegrationPointCoordinates(xyz, distances);
232 for (
int ip = 0; ip < static_cast<int>(distances.size()); ++ip)
239 _ip_data[k].non_local_assemblers.push_back(
240 {la->getIPDataPtr(ip),
241 std::numeric_limits<double>::quiet_NaN(),
245 if (
_ip_data[k].non_local_assemblers.size() == 0)
250 double a_k_sum_m = 0;
251 for (
auto const& tuple :
_ip_data[k].non_local_assemblers)
253 double const distance2_m = tuple.distance2;
255 auto const& w_m = tuple.ip_l_pointer->integration_weight;
257 a_k_sum_m += w_m *
alpha_0(distance2_m);
264 for (
auto& tuple :
_ip_data[k].non_local_assemblers)
266 double const distance2_l = tuple.distance2;
267 double const a_kl =
alpha_0(distance2_l) / a_k_sum_m;
271 auto const w_l = tuple.ip_l_pointer->integration_weight;
272 tuple.alpha_kl_times_w_l = a_kl * w_l;
278 int integration_point)
const
282 Eigen::Vector3d xyz = Eigen::Vector3d::Zero();
284 for (
int i = 0; i < N.size(); ++i)
286 auto const& node_coordinates{nodes[i]->asEigenVector3d()};
287 xyz += node_coordinates * N[i];
296 Eigen::Vector3d
const& coords,
297 std::vector<double>& distances)
const override
299 unsigned const n_integration_points =
302 distances.resize(n_integration_points);
304 for (
unsigned ip = 0; ip < n_integration_points; ip++)
306 auto const& xyz =
_ip_data[ip].coordinates;
307 distances[ip] = (xyz - coords).squaredNorm();
312 std::vector<double>
const& ,
313 std::vector<double>
const& ,
314 std::vector<double>& ,
315 std::vector<double>& ,
316 std::vector<double>& )
override
319 "SmallDeformationNonlocalLocalAssembler: assembly without jacobian "
325 std::vector<double>
const& local_x)
override
327 auto const n_integration_points =
335 for (
unsigned ip = 0; ip < n_integration_points; ip++)
338 auto const& dNdx =
_ip_data[ip].dNdx;
344 DisplacementDim, ShapeFunction::NPOINTS,
347 auto const& eps_prev =
_ip_data[ip].eps_prev;
348 auto const& sigma_prev =
_ip_data[ip].sigma_prev;
353 auto& state =
_ip_data[ip].material_state_variables;
354 double const& damage_prev =
_ip_data[ip].damage_prev;
358 Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
359 local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
366 DisplacementDim>::MaterialStateVariables>
371 KelvinVectorType
const sigma_eff_prev =
376 variables_prev.
stress.emplace<
387 auto&& solution =
_ip_data[ip].solid_material.integrateStress(
388 variables_prev, variables, t, x_position, dt, *state);
392 OGS_FATAL(
"Computation of local constitutive relation failed.");
395 std::tie(sigma, state, C) = std::move(*solution);
399 auto const& ehlers_material =
401 DisplacementDim
> const&>(
_ip_data[ip].solid_material);
402 auto const damage_properties =
404 auto const material_properties =
405 ehlers_material.evaluatedMaterialProperties(t, x_position);
411 *
_ip_data[ip].material_state_variables);
413 double const eps_p_eff_diff =
414 state_vars.
eps_p.eff - state_vars.eps_p_prev.eff;
417 eps_p_eff_diff, sigma,
_ip_data[ip].kappa_d_prev,
418 damage_properties.h_d, material_properties);
425 for (
auto const& tuple :
429 tuple.ip_l_pointer->activated =
true;
438 std::vector<double>
const& local_x,
439 std::vector<double>
const& ,
440 std::vector<double>& local_b_data,
441 std::vector<double>& local_Jac_data)
override
443 auto const local_matrix_size = local_x.size();
446 local_Jac_data, local_matrix_size, local_matrix_size);
449 local_b_data, local_matrix_size);
451 unsigned const n_integration_points =
458 for (
unsigned ip = 0; ip < n_integration_points; ip++)
460 auto const& w =
_ip_data[ip].integration_weight;
463 auto const& dNdx =
_ip_data[ip].dNdx;
469 DisplacementDim, ShapeFunction::NPOINTS,
475 double& damage =
_ip_data[ip].damage;
478 double nonlocal_kappa_d = 0;
482 for (
auto const& tuple :
_ip_data[ip].non_local_assemblers)
485 double const kappa_d_l = tuple.ip_l_pointer->kappa_d;
486 double const a_kl_times_w_l = tuple.alpha_kl_times_w_l;
487 nonlocal_kappa_d += a_kl_times_w_l * kappa_d_l;
491 auto const& ehlers_material =
493 DisplacementDim
> const&>(
_ip_data[ip].solid_material);
508 nonlocal_kappa_d = std::max(0., nonlocal_kappa_d);
512 auto const damage_properties =
513 ehlers_material.evaluatedDamageProperties(t,
516 damage_properties.alpha_d,
517 damage_properties.beta_d);
518 damage = std::max(0., damage);
520 sigma = sigma * (1. - damage);
523 local_b.noalias() -= B.transpose() * sigma * w;
524 local_Jac.noalias() += B.transpose() * C * (1. - damage) * B * w;
530 unsigned const n_integration_points =
533 for (
unsigned ip = 0; ip < n_integration_points; ip++)
540 Eigen::VectorXd
const& ,
541 double const ,
double const ,
544 unsigned const n_integration_points =
547 for (
unsigned ip = 0; ip < n_integration_points; ip++)
556 double& crack_volume)
override
559 auto local_x = x.
get(indices);
561 auto u = Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
562 local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
564 int const n_integration_points =
567 for (
int ip = 0; ip < n_integration_points; ip++)
569 auto const& dNdx =
_ip_data[ip].dNdx;
570 auto const& d =
_ip_data[ip].damage;
571 auto const& w =
_ip_data[ip].integration_weight;
574 (dNdx.template topRows<DisplacementDim>() *
575 u.reshaped(ShapeFunction::NPOINTS, DisplacementDim))
577 crack_volume += div_u * d * w;
582 const unsigned integration_point)
const override
587 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
591 std::vector<double>& nodal_values)
const override
593 nodal_values.clear();
595 nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
597 unsigned const n_integration_points =
600 for (
unsigned ip = 0; ip < n_integration_points; ip++)
602 auto const& w =
_ip_data[ip].integration_weight;
605 auto const& dNdx =
_ip_data[ip].dNdx;
611 DisplacementDim, ShapeFunction::NPOINTS,
616 local_b.noalias() += B.transpose() * sigma * w;
624 std::vector<GlobalVector*>
const& ,
625 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
626 std::vector<double>& cache)
const override
632 [](
auto const& ip_data)
633 {
return ip_data.free_energy_density; });
640 std::vector<GlobalVector*>
const& ,
641 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
642 std::vector<double>& cache)
const override
648 [](
auto const& ip_data) {
return *ip_data.eps_p_V; });
655 std::vector<GlobalVector*>
const& ,
656 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
657 std::vector<double>& cache)
const override
663 [](
auto const& ip_data) {
return *ip_data.eps_p_D_xx; });
670 std::vector<GlobalVector*>
const& ,
671 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
672 std::vector<double>& cache)
const override
680 std::vector<GlobalVector*>
const& ,
681 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
682 std::vector<double>& cache)
const override
707 ip_data.kappa_d = value;
712 unsigned const n_integration_points =
715 std::vector<double> result_values;
716 result_values.resize(n_integration_points);
718 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
720 result_values[ip] =
_ip_data[ip].kappa_d;
723 return result_values;
728 std::vector<GlobalVector*>
const& ,
729 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
730 std::vector<double>& cache)
const override
742 DisplacementDim>::MaterialStateVariables
const&
745 return *
_ip_data[integration_point].material_state_variables;
750 std::size_t
const component)
const
755 for (
auto const& ip_data :
_ip_data)
759 cache.push_back(ip_data.sigma[component]);
763 cache.push_back(ip_data.sigma[component] / std::sqrt(2));
771 std::vector<double>& cache, std::size_t
const component)
const
776 for (
auto const& ip_data :
_ip_data)
779 cache.push_back(ip_data.eps[component]);
781 cache.push_back(ip_data.eps[component] / std::sqrt(2));
795 std::vector<IpData, Eigen::aligned_allocator<IpData>>
_ip_data;
803 ShapeFunction::NPOINTS * DisplacementDim;
DamageProperties evaluatedDamageProperties(double const t, ParameterLib::SpatialPosition const &x) const
KelvinVector mechanical_strain
Global vector based on Eigen vector.
double get(IndexType rowId) const
get entry
constexpr double getWeight() const
virtual Node *const * getNodes() const =0
Get array of element nodes.
std::size_t getID() const
Returns the ID of the element.
unsigned getIntegrationOrder() const
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
unsigned getNumberOfPoints() const
void setElementID(std::size_t element_id)
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
MatrixType< _number_of_dof, _number_of_dof > StiffnessMatrixType
VectorType< _number_of_dof > NodalForceVectorType
Rhs residual.
VectorType< _kelvin_vector_size > KelvinVectorType
auto & selectSolidConstitutiveRelation(SolidMaterialsMap const &constitutive_relations, MeshLib::PropertyVector< int > const *const material_ids, std::size_t const element_id)
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), kelvin_vector_dimensions(DisplacementDim), Eigen::RowMajor > KelvinMatrixType
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
std::vector< std::size_t > findElementsWithinRadius(Element const &start_element, double const radius_squared)
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
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)
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.
std::vector< double > const & getIntegrationPointScalarData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
std::vector< double > const & getIntegrationPointKelvinVectorData(IntegrationPointDataVector const &ip_data_vector, MemberType IpData::*const member, std::vector< double > &cache)
std::size_t setIntegrationPointKelvinVectorData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)
std::size_t setIntegrationPointScalarData(double const *values, IntegrationPointDataVector &ip_data_vector, MemberType IpData::*const member)
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
VectorType< GlobalDim > GlobalDimVectorType
VectorType< ShapeFunction::NPOINTS > NodalVectorType
PlasticStrain< KelvinVector > eps_p
plastic part of the state.