39 namespace SmallDeformationNonlocal
45 template <
typename ShapeMatrixType>
48 std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>>
N;
51 template <
typename ShapeFunction,
typename IntegrationMethod,
81 bool const is_axially_symmetric,
82 unsigned const integration_order,
89 unsigned const n_integration_points =
92 _ip_data.reserve(n_integration_points);
95 auto const shape_matrices =
97 DisplacementDim>(e, is_axially_symmetric,
100 auto& solid_material =
105 auto* ehlers_solid_material =
dynamic_cast<
108 if (ehlers_solid_material ==
nullptr)
111 "The SmallDeformationNonlocalLocal process supports only "
112 "Ehlers material at the moment. For other materials the "
113 "interface must be extended first.");
116 for (
unsigned ip = 0; ip < n_integration_points; ip++)
118 _ip_data.emplace_back(*ehlers_solid_material);
120 auto const& sm = shape_matrices[ip];
123 sm.integralMeasure * sm.detJ;
126 ip_data.dNdx = sm.dNdx;
129 ip_data.sigma.setZero(
136 ip_data.sigma_prev.resize(
139 ip_data.eps_prev.resize(
150 double const* values,
151 int const integration_order)
override
153 if (integration_order !=
157 "Setting integration point initial conditions; The integration "
158 "order of the local assembler for element {:d} is different "
159 "from the integration order in the initial condition.",
163 if (
name ==
"sigma_ip")
168 if (
name ==
"kappa_d_ip")
178 std::string
const&
name, std::vector<double>
const& value)
override
180 if (
name ==
"kappa_d_ip")
182 if (value.size() != 1)
185 "CellData for kappa_d initial conditions has wrong number "
186 "of components. 1 expected, got {:d}.",
195 double const internal_length2 =
_process_data.internal_length_squared;
196 return (distance2 > internal_length2)
198 : (1 - distance2 / (internal_length2)) *
199 (1 - distance2 / (internal_length2));
206 DisplacementDim>>>
const& local_assemblers)
override
211 unsigned const n_integration_points =
214 std::vector<double> distances;
220 for (
unsigned k = 0; k < n_integration_points; k++)
226 auto const& xyz =
_ip_data[k].coordinates;
229 for (
auto const search_element_id : search_element_ids)
231 auto const& la = local_assemblers[search_element_id];
232 la->getIntegrationPointCoordinates(xyz, distances);
233 for (
int ip = 0; ip < static_cast<int>(distances.size()); ++ip)
240 _ip_data[k].non_local_assemblers.push_back(
241 {la->getIPDataPtr(ip),
242 std::numeric_limits<double>::quiet_NaN(),
246 if (
_ip_data[k].non_local_assemblers.size() == 0)
251 double a_k_sum_m = 0;
252 for (
auto const& tuple :
_ip_data[k].non_local_assemblers)
254 double const distance2_m = tuple.distance2;
256 auto const& w_m = tuple.ip_l_pointer->integration_weight;
258 a_k_sum_m += w_m *
alpha_0(distance2_m);
265 for (
auto& tuple :
_ip_data[k].non_local_assemblers)
267 double const distance2_l = tuple.distance2;
268 double const a_kl =
alpha_0(distance2_l) / a_k_sum_m;
272 auto const w_l = tuple.ip_l_pointer->integration_weight;
273 tuple.alpha_kl_times_w_l = a_kl * w_l;
279 int integration_point)
const
283 Eigen::Vector3d xyz = Eigen::Vector3d::Zero();
285 for (
int i = 0; i < N.size(); ++i)
287 Eigen::Map<Eigen::Vector3d const> node_coordinates{
289 xyz += node_coordinates * N[i];
298 Eigen::Vector3d
const& coords,
299 std::vector<double>& distances)
const override
301 unsigned const n_integration_points =
304 distances.resize(n_integration_points);
306 for (
unsigned ip = 0; ip < n_integration_points; ip++)
308 auto const& xyz =
_ip_data[ip].coordinates;
309 distances[ip] = (xyz - coords).squaredNorm();
314 std::vector<double>
const& ,
315 std::vector<double>
const& ,
316 std::vector<double>& ,
317 std::vector<double>& ,
318 std::vector<double>& )
override
321 "SmallDeformationNonlocalLocalAssembler: assembly without jacobian "
327 std::vector<double>
const& local_x)
override
329 auto const n_integration_points =
337 for (
unsigned ip = 0; ip < n_integration_points; ip++)
342 auto const& dNdx =
_ip_data[ip].dNdx;
348 DisplacementDim, ShapeFunction::NPOINTS,
351 auto const& eps_prev =
_ip_data[ip].eps_prev;
352 auto const& sigma_prev =
_ip_data[ip].sigma_prev;
357 auto& state =
_ip_data[ip].material_state_variables;
358 double const& damage_prev =
_ip_data[ip].damage_prev;
362 Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
363 local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
370 DisplacementDim>::MaterialStateVariables>
380 variables_prev[
static_cast<int>(MPL::Variable::stress)]
384 variables_prev[
static_cast<int>(MPL::Variable::mechanical_strain)]
390 variables[
static_cast<int>(MPL::Variable::mechanical_strain)]
397 auto&& solution =
_ip_data[ip].solid_material.integrateStress(
398 variables_prev, variables, t, x_position, dt, *state);
402 OGS_FATAL(
"Computation of local constitutive relation failed.");
405 std::tie(sigma, state, C) = std::move(*solution);
409 auto const& ehlers_material =
411 DisplacementDim
> const&>(
_ip_data[ip].solid_material);
412 auto const damage_properties =
414 auto const material_properties =
415 ehlers_material.evaluatedMaterialProperties(t, x_position);
421 *
_ip_data[ip].material_state_variables);
423 double const eps_p_eff_diff =
424 state_vars.
eps_p.eff - state_vars.eps_p_prev.eff;
426 _ip_data[ip].kappa_d = calculateDamageKappaD<DisplacementDim>(
427 eps_p_eff_diff, sigma,
_ip_data[ip].kappa_d_prev,
428 damage_properties.h_d, material_properties);
435 for (
auto const& tuple :
439 tuple.ip_l_pointer->activated =
true;
448 std::vector<double>
const& local_x,
449 std::vector<double>
const& ,
450 const double ,
const double ,
451 std::vector<double>& ,
452 std::vector<double>& ,
453 std::vector<double>& local_b_data,
454 std::vector<double>& local_Jac_data)
override
456 auto const local_matrix_size = local_x.size();
458 auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
459 local_Jac_data, local_matrix_size, local_matrix_size);
461 auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
462 local_b_data, local_matrix_size);
464 unsigned const n_integration_points =
471 for (
unsigned ip = 0; ip < n_integration_points; ip++)
474 auto const& w =
_ip_data[ip].integration_weight;
477 auto const& dNdx =
_ip_data[ip].dNdx;
483 DisplacementDim, ShapeFunction::NPOINTS,
489 double& damage =
_ip_data[ip].damage;
492 double nonlocal_kappa_d = 0;
496 for (
auto const& tuple :
_ip_data[ip].non_local_assemblers)
499 double const kappa_d_l = tuple.ip_l_pointer->kappa_d;
500 double const a_kl_times_w_l = tuple.alpha_kl_times_w_l;
501 nonlocal_kappa_d += a_kl_times_w_l * kappa_d_l;
505 auto const& ehlers_material =
507 DisplacementDim
> const&>(
_ip_data[ip].solid_material);
522 nonlocal_kappa_d = std::max(0., nonlocal_kappa_d);
526 auto const damage_properties =
527 ehlers_material.evaluatedDamageProperties(t,
530 damage_properties.alpha_d,
531 damage_properties.beta_d);
532 damage = std::max(0., damage);
534 sigma = sigma * (1. - damage);
537 local_b.noalias() -= B.transpose() * sigma * w;
538 local_Jac.noalias() += B.transpose() * C * (1. - damage) * B * w;
544 unsigned const n_integration_points =
547 for (
unsigned ip = 0; ip < n_integration_points; ip++)
554 double const ,
double const )
override
556 unsigned const n_integration_points =
559 for (
unsigned ip = 0; ip < n_integration_points; ip++)
568 double& crack_volume)
override
571 auto local_x = x.
get(indices);
573 auto u = Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
574 local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
576 int const n_integration_points =
579 for (
int ip = 0; ip < n_integration_points; ip++)
581 auto const& dNdx =
_ip_data[ip].dNdx;
582 auto const& d =
_ip_data[ip].damage;
583 auto const& w =
_ip_data[ip].integration_weight;
587 ShapeFunction::NPOINTS>(u, dNdx);
588 crack_volume += div_u * d * w;
593 const unsigned integration_point)
const override
598 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
602 std::vector<double>& nodal_values)
const override
604 nodal_values.clear();
605 auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
606 nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
608 unsigned const n_integration_points =
611 for (
unsigned ip = 0; ip < n_integration_points; ip++)
613 auto const& w =
_ip_data[ip].integration_weight;
616 auto const& dNdx =
_ip_data[ip].dNdx;
622 DisplacementDim, ShapeFunction::NPOINTS,
627 local_b.noalias() += B.transpose() * sigma * w;
635 std::vector<GlobalVector*>
const& ,
636 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
637 std::vector<double>& cache)
const override
644 [](
auto const& ip_data) {
return ip_data.free_energy_density; });
651 std::vector<GlobalVector*>
const& ,
652 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
653 std::vector<double>& cache)
const override
659 [](
auto const& ip_data) {
return *ip_data.eps_p_V; });
666 std::vector<GlobalVector*>
const& ,
667 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
668 std::vector<double>& cache)
const override
674 [](
auto const& ip_data) {
return *ip_data.eps_p_D_xx; });
681 std::vector<GlobalVector*>
const& ,
682 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
683 std::vector<double>& cache)
const override
685 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
691 std::vector<GlobalVector*>
const& ,
692 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
693 std::vector<double>& cache)
const override
695 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
701 return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
710 return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
718 ip_data.kappa_d = value;
723 unsigned const n_integration_points =
726 std::vector<double> result_values;
727 result_values.resize(n_integration_points);
729 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
731 result_values[ip] =
_ip_data[ip].kappa_d;
734 return result_values;
739 std::vector<GlobalVector*>
const& ,
740 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& ,
741 std::vector<double>& cache)
const override
753 DisplacementDim>::MaterialStateVariables
const&
756 return *
_ip_data[integration_point].material_state_variables;
761 std::size_t
const component)
const
766 for (
auto const& ip_data :
_ip_data)
770 cache.push_back(ip_data.sigma[component]);
774 cache.push_back(ip_data.sigma[component] / std::sqrt(2));
782 std::vector<double>& cache, std::size_t
const component)
const
787 for (
auto const& ip_data :
_ip_data)
790 cache.push_back(ip_data.eps[component]);
792 cache.push_back(ip_data.eps[component] / std::sqrt(2));
807 std::vector<IpData, Eigen::aligned_allocator<IpData>>
_ip_data;
815 ShapeFunction::NPOINTS * DisplacementDim;
DamageProperties evaluatedDamageProperties(double const t, ParameterLib::SpatialPosition const &x) const
Global vector based on Eigen vector.
double get(IndexType rowId) const
get entry
const T * getCoords() const
virtual Node *const * getNodes() const =0
Get array of element nodes.
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
MatrixType< _number_of_dof, _number_of_dof > StiffnessMatrixType
VectorType< _number_of_dof > NodalForceVectorType
Rhs residual.
VectorType< _kelvin_vector_size > KelvinVectorType
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< 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(std::vector< IntegrationPointData, Eigen::aligned_allocator< IntegrationPointData >> const &ip_data, MemberType member, std::vector< double > &cache)
std::size_t setIntegrationPointScalarData(double const *values, std::vector< IntegrationPointData, Eigen::aligned_allocator< IntegrationPointData >> &ip_data, MemberType member)
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
VectorType< GlobalDim > GlobalDimVectorType
VectorType< ShapeFunction::NPOINTS > NodalVectorType
PlasticStrain< KelvinVector > eps_p
plastic part of the state.
std::vector< ShapeMatrixType, Eigen::aligned_allocator< ShapeMatrixType > > N