40 namespace SmallDeformationNonlocal
44 template <
typename ShapeMatrixType>
47 std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>>
N;
50 template <
typename ShapeFunction,
typename IntegrationMethod,
78 bool const is_axially_symmetric,
79 unsigned const integration_order,
81 : _process_data(process_data),
82 _integration_method(integration_order),
84 _is_axially_symmetric(is_axially_symmetric)
86 unsigned const n_integration_points =
87 _integration_method.getNumberOfPoints();
89 _ip_data.reserve(n_integration_points);
90 _secondary_data.N.resize(n_integration_points);
92 auto const shape_matrices =
94 IntegrationMethod, DisplacementDim>(
95 e, is_axially_symmetric, _integration_method);
97 auto& solid_material =
99 _process_data.solid_materials,
100 _process_data.material_ids,
102 auto* ehlers_solid_material =
dynamic_cast< 105 if (ehlers_solid_material ==
nullptr)
108 "The SmallDeformationNonlocalLocal process supports only " 109 "Ehlers material at the moment. For other materials the " 110 "interface must be extended first.");
113 for (
unsigned ip = 0; ip < n_integration_points; ip++)
115 _ip_data.emplace_back(*ehlers_solid_material);
116 auto& ip_data = _ip_data[ip];
117 auto const& sm = shape_matrices[ip];
118 _ip_data[ip].integration_weight =
119 _integration_method.getWeightedPoint(ip).getWeight() *
120 sm.integralMeasure * sm.detJ;
123 ip_data.dNdx = sm.dNdx;
127 DisplacementDim>::value);
129 DisplacementDim>::value);
132 ip_data.sigma_prev.resize(
134 DisplacementDim>::value);
135 ip_data.eps_prev.resize(
137 DisplacementDim>::value);
139 _secondary_data.N[ip] = shape_matrices[ip].N;
141 ip_data.coordinates = getSingleIntegrationPointCoordinates(ip);
146 double const* values,
147 int const integration_order)
override 149 if (integration_order !=
150 static_cast<int>(_integration_method.getIntegrationOrder()))
153 "Setting integration point initial conditions; The integration " 154 "order of the local assembler for element %d is different from " 155 "the integration order in the initial condition.",
159 if (name ==
"sigma_ip")
161 return setSigma(values);
164 if (name ==
"kappa_d_ip")
166 return setKappaD(values);
173 std::string
const&
name, std::vector<double>
const& value)
override 175 if (name ==
"kappa_d_ip")
177 if (value.size() != 1)
180 "CellData for kappa_d initial conditions has wrong number " 181 "of components. 1 expected, got %d.",
190 double const internal_length2 = _process_data.internal_length_squared;
191 return (distance2 > internal_length2)
193 : (1 - distance2 / (internal_length2)) *
194 (1 - distance2 / (internal_length2));
201 DisplacementDim>>>
const& local_assemblers)
override 204 _element, _process_data.internal_length_squared);
206 unsigned const n_integration_points =
207 _integration_method.getNumberOfPoints();
209 std::vector<double> distances;
215 for (
unsigned k = 0; k < n_integration_points; k++)
221 auto const& xyz = _ip_data[k].coordinates;
224 for (
auto const search_element_id : search_element_ids)
226 auto const& la = local_assemblers[search_element_id];
227 la->getIntegrationPointCoordinates(xyz, distances);
228 for (
int ip = 0; ip < static_cast<int>(distances.size()); ++ip)
230 if (distances[ip] >= _process_data.internal_length_squared)
233 _ip_data[k].non_local_assemblers.push_back(
234 {la->getIPDataPtr(ip),
235 std::numeric_limits<double>::quiet_NaN(),
239 if (_ip_data[k].non_local_assemblers.size() == 0)
244 double a_k_sum_m = 0;
245 for (
auto const& tuple : _ip_data[k].non_local_assemblers)
247 double const distance2_m = tuple.distance2;
249 auto const& w_m = tuple.ip_l_pointer->integration_weight;
251 a_k_sum_m += w_m * alpha_0(distance2_m);
258 for (
auto& tuple : _ip_data[k].non_local_assemblers)
260 double const distance2_l = tuple.distance2;
261 double const a_kl = alpha_0(distance2_l) / a_k_sum_m;
265 auto const w_l = tuple.ip_l_pointer->integration_weight;
266 tuple.alpha_kl_times_w_l = a_kl * w_l;
272 int integration_point)
const 274 auto const&
N = _secondary_data.N[integration_point];
276 Eigen::Vector3d xyz = Eigen::Vector3d::Zero();
277 auto* nodes = _element.getNodes();
278 for (
int i = 0; i <
N.size(); ++i)
280 Eigen::Map<Eigen::Vector3d const> node_coordinates{
281 nodes[i]->getCoords(), 3};
282 xyz += node_coordinates *
N[i];
291 Eigen::Vector3d
const& coords,
292 std::vector<double>& distances)
const override 294 unsigned const n_integration_points =
295 _integration_method.getNumberOfPoints();
297 distances.resize(n_integration_points);
299 for (
unsigned ip = 0; ip < n_integration_points; ip++)
301 auto const& xyz = _ip_data[ip].coordinates;
302 distances[ip] = (xyz - coords).squaredNorm();
306 void assemble(
double const , std::vector<double>
const& ,
307 std::vector<double>& ,
308 std::vector<double>& ,
309 std::vector<double>& )
override 312 "SmallDeformationNonlocalLocalAssembler: assembly without jacobian " 318 std::vector<double>
const& local_x)
override 320 auto const n_integration_points =
321 _integration_method.getNumberOfPoints();
326 for (
unsigned ip = 0; ip < n_integration_points; ip++)
330 auto const&
N = _ip_data[ip].N;
331 auto const& dNdx = _ip_data[ip].dNdx;
334 interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
337 DisplacementDim, ShapeFunction::NPOINTS,
339 _is_axially_symmetric);
340 auto const& eps_prev = _ip_data[ip].eps_prev;
341 auto const& sigma_prev = _ip_data[ip].sigma_prev;
343 auto& eps = _ip_data[ip].eps;
344 auto& sigma = _ip_data[ip].sigma;
345 auto& C = _ip_data[ip].C;
346 auto& state = _ip_data[ip].material_state_variables;
347 double const& damage_prev = _ip_data[ip].damage_prev;
351 Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
352 local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
359 DisplacementDim>::MaterialStateVariables>
369 auto&& solution = _ip_data[ip].solid_material.integrateStress(
370 t, x_position, _process_data.dt, eps_prev, eps, sigma_eff_prev,
371 *state, _process_data.reference_temperature);
374 OGS_FATAL(
"Computation of local constitutive relation failed.");
376 std::tie(sigma, state, C) = std::move(*solution);
380 auto const& ehlers_material =
382 DisplacementDim
> const&>(_ip_data[ip].solid_material);
383 auto const damage_properties =
385 auto const material_properties =
386 ehlers_material.evaluatedMaterialProperties(t, x_position);
392 *_ip_data[ip].material_state_variables);
394 double const eps_p_eff_diff =
395 state_vars.
eps_p.eff - state_vars.eps_p_prev.eff;
397 _ip_data[ip].kappa_d = calculateDamageKappaD<DisplacementDim>(
398 eps_p_eff_diff, sigma, _ip_data[ip].kappa_d_prev,
399 damage_properties.h_d, material_properties);
401 if (!_ip_data[ip].active_self)
403 _ip_data[ip].active_self |= _ip_data[ip].kappa_d > 0;
404 if (_ip_data[ip].active_self)
406 for (
auto const& tuple :
407 _ip_data[ip].non_local_assemblers)
410 tuple.ip_l_pointer->activated =
true;
419 std::vector<double>
const& local_x,
420 std::vector<double>
const& ,
421 const double ,
const double ,
422 std::vector<double>& ,
423 std::vector<double>& ,
424 std::vector<double>& local_b_data,
425 std::vector<double>& local_Jac_data)
override 427 auto const local_matrix_size = local_x.size();
429 auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
430 local_Jac_data, local_matrix_size, local_matrix_size);
432 auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
433 local_b_data, local_matrix_size);
435 unsigned const n_integration_points =
436 _integration_method.getNumberOfPoints();
442 for (
unsigned ip = 0; ip < n_integration_points; ip++)
445 auto const& w = _ip_data[ip].integration_weight;
447 auto const&
N = _ip_data[ip].N;
448 auto const& dNdx = _ip_data[ip].dNdx;
451 interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
454 DisplacementDim, ShapeFunction::NPOINTS,
456 _is_axially_symmetric);
458 auto& sigma = _ip_data[ip].sigma;
459 auto& C = _ip_data[ip].C;
460 double& damage = _ip_data[ip].damage;
463 double nonlocal_kappa_d = 0;
465 if (_ip_data[ip].active_self || _ip_data[ip].activated)
467 for (
auto const& tuple : _ip_data[ip].non_local_assemblers)
470 double const kappa_d_l = tuple.ip_l_pointer->kappa_d;
471 double const a_kl_times_w_l = tuple.alpha_kl_times_w_l;
472 nonlocal_kappa_d += a_kl_times_w_l * kappa_d_l;
476 auto const& ehlers_material =
478 DisplacementDim
> const&>(_ip_data[ip].solid_material);
493 nonlocal_kappa_d = std::max(0., nonlocal_kappa_d);
497 auto const damage_properties =
498 ehlers_material.evaluatedDamageProperties(t,
501 damage_properties.alpha_d,
502 damage_properties.beta_d);
503 damage = std::max(0., damage);
505 sigma = sigma * (1. - damage);
508 local_b.noalias() -= B.transpose() * sigma * w;
509 local_Jac.noalias() += B.transpose() * C * (1. - damage) * B * w;
515 double const )
override 517 unsigned const n_integration_points =
518 _integration_method.getNumberOfPoints();
520 for (
unsigned ip = 0; ip < n_integration_points; ip++)
522 _ip_data[ip].pushBackState();
528 GlobalVector
const& x,
529 double& crack_volume)
override 532 auto local_x = x.get(indices);
534 auto u = Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
535 local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
537 int const n_integration_points =
538 _integration_method.getNumberOfPoints();
540 for (
int ip = 0; ip < n_integration_points; ip++)
542 auto const& dNdx = _ip_data[ip].dNdx;
543 auto const& d = _ip_data[ip].damage;
544 auto const& w = _ip_data[ip].integration_weight;
548 ShapeFunction::NPOINTS>(
u, dNdx);
549 crack_volume += div_u * d * w;
554 const unsigned integration_point)
const override 556 auto const&
N = _secondary_data.N[integration_point];
559 return Eigen::Map<const Eigen::RowVectorXd>(
N.data(),
N.size());
563 std::vector<double>& nodal_values)
const override 565 nodal_values.clear();
566 auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
567 nodal_values, ShapeFunction::NPOINTS * DisplacementDim);
569 unsigned const n_integration_points =
570 _integration_method.getNumberOfPoints();
572 for (
unsigned ip = 0; ip < n_integration_points; ip++)
574 auto const& w = _ip_data[ip].integration_weight;
576 auto const&
N = _ip_data[ip].N;
577 auto const& dNdx = _ip_data[ip].dNdx;
580 interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
583 DisplacementDim, ShapeFunction::NPOINTS,
585 _is_axially_symmetric);
586 auto& sigma = _ip_data[ip].sigma;
588 local_b.noalias() += B.transpose() * sigma * w;
596 GlobalVector
const& ,
598 std::vector<double>& cache)
const override 601 cache.reserve(_ip_data.size());
603 for (
auto const& ip_data : _ip_data)
605 cache.push_back(ip_data.free_energy_density);
613 GlobalVector
const& ,
615 std::vector<double>& cache)
const override 618 cache.reserve(_ip_data.size());
620 for (
auto const& ip_data : _ip_data)
622 cache.push_back(*ip_data.eps_p_V);
630 GlobalVector
const& ,
632 std::vector<double>& cache)
const override 635 cache.reserve(_ip_data.size());
637 for (
auto const& ip_data : _ip_data)
639 cache.push_back(*ip_data.eps_p_D_xx);
647 GlobalVector
const& ,
649 std::vector<double>& cache)
const override 651 static const int kelvin_vector_size =
653 DisplacementDim>::value;
654 auto const num_intpts = _ip_data.size();
658 double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
659 cache, kelvin_vector_size, num_intpts);
661 for (
unsigned ip = 0; ip < num_intpts; ++ip)
663 auto const& sigma = _ip_data[ip].sigma;
673 GlobalVector
const& ,
675 std::vector<double>& cache)
const override 677 auto const kelvin_vector_size =
679 DisplacementDim>::value;
680 auto const num_intpts = _ip_data.size();
684 double, kelvin_vector_size, Eigen::Dynamic, Eigen::RowMajor>>(
685 cache, kelvin_vector_size, num_intpts);
687 for (
unsigned ip = 0; ip < num_intpts; ++ip)
689 auto const& eps = _ip_data[ip].eps;
699 auto const kelvin_vector_size =
701 DisplacementDim>::value;
702 auto const n_integration_points = _ip_data.size();
704 std::vector<double> ip_sigma_values;
706 Eigen::Map<Eigen::Matrix<double, kelvin_vector_size, Eigen::Dynamic,
707 Eigen::ColMajor>
const>(
708 values, kelvin_vector_size, n_integration_points);
710 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
714 sigma_values.col(ip));
717 return n_integration_points;
725 auto const kelvin_vector_size =
727 DisplacementDim>::value;
728 auto const n_integration_points = _ip_data.size();
730 std::vector<double> ip_sigma_values;
732 double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
733 ip_sigma_values, n_integration_points, kelvin_vector_size);
735 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
737 auto const& sigma = _ip_data[ip].sigma;
742 return ip_sigma_values;
747 unsigned const n_integration_points =
748 _integration_method.getNumberOfPoints();
750 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
752 _ip_data[ip].kappa_d = values[ip];
754 return n_integration_points;
759 for (
auto& ip_data : _ip_data)
761 ip_data.kappa_d = value;
766 unsigned const n_integration_points =
767 _integration_method.getNumberOfPoints();
769 std::vector<double> result_values;
770 result_values.resize(n_integration_points);
772 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
774 result_values[ip] = _ip_data[ip].kappa_d;
777 return result_values;
782 GlobalVector
const& ,
784 std::vector<double>& cache)
const override 787 cache.reserve(_ip_data.size());
789 for (
auto const& ip_data : _ip_data)
791 cache.push_back(ip_data.damage);
799 return _integration_method.getNumberOfPoints();
803 DisplacementDim>::MaterialStateVariables
const&
806 return *_ip_data[integration_point].material_state_variables;
811 std::size_t
const component)
const 814 cache.reserve(_ip_data.size());
816 for (
auto const& ip_data : _ip_data)
819 cache.push_back(ip_data.sigma[component]);
821 cache.push_back(ip_data.sigma[component] / std::sqrt(2));
828 std::vector<double>& cache, std::size_t
const component)
const 831 cache.reserve(_ip_data.size());
833 for (
auto const& ip_data : _ip_data)
836 cache.push_back(ip_data.eps[component]);
838 cache.push_back(ip_data.eps[component] / std::sqrt(2));
847 return &_ip_data[ip];
864 static const int displacement_size =
865 ShapeFunction::NPOINTS * DisplacementDim;
Eigen::Matrix< double, Eigen::MatrixBase< Derived >::RowsAtCompileTime, 1 > symmetricTensorToKelvinVector(Eigen::MatrixBase< Derived > const &v)
Kelvin vector dimensions for given displacement dimension.
MatrixType< _number_of_dof, _number_of_dof > StiffnessMatrixType
VectorType< _number_of_dof > NodalForceVectorType
Rhs residual.
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool is_axially_symmetric, IntegrationMethod const &integration_method)
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)
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, KelvinVectorDimensions< DisplacementDim >::value, Eigen::RowMajor > KelvinMatrixType
void setElementID(std::size_t element_id)
PlasticStrain< KelvinVector > eps_p
plastic part of the state.
DamageProperties evaluatedDamageProperties(double const t, ProcessLib::SpatialPosition const &x) const
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
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.
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
VectorType< ShapeFunction::NPOINTS > NodalVectorType
VectorType< GlobalDim > GlobalDimVectorType
std::vector< std::size_t > findElementsWithinRadius(Element const &start_element, double const radius_squared)
#define OGS_FATAL(fmt,...)
void setIntegrationPoint(unsigned integration_point)
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, 1, Eigen::ColMajor > KelvinVectorType