OGS
ThermoRichardsMechanicsFEM-impl.h
Go to the documentation of this file.
1
12#pragma once
13
14#include <spdlog/fmt/bundled/format.h>
15
16#include <Eigen/LU>
17#include <cassert>
18
29
30namespace ProcessLib
31{
32namespace ThermoRichardsMechanics
33{
34template <typename ShapeFunctionDisplacement, typename ShapeFunction,
35 int DisplacementDim, typename ConstitutiveTraits>
36ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunction,
37 DisplacementDim, ConstitutiveTraits>::
38 ThermoRichardsMechanicsLocalAssembler(
39 MeshLib::Element const& e,
40 std::size_t const /*local_matrix_size*/,
41 NumLib::GenericIntegrationMethod const& integration_method,
42 bool const is_axially_symmetric,
44 process_data)
45 : LocalAssemblerInterface<DisplacementDim, ConstitutiveTraits>(
46 e, integration_method, is_axially_symmetric, process_data)
47{
48 unsigned const n_integration_points =
50
51 ip_data_.resize(n_integration_points);
52
53 auto const shape_matrices_u =
54 NumLib::initShapeMatrices<ShapeFunctionDisplacement,
56 DisplacementDim>(e, is_axially_symmetric,
57 integration_method);
58
59 auto const shape_matrices =
61 DisplacementDim>(e, is_axially_symmetric,
62 integration_method);
63
64 for (unsigned ip = 0; ip < n_integration_points; ip++)
65 {
66 auto& ip_data = ip_data_[ip];
67 auto const& sm_u = shape_matrices_u[ip];
68 ip_data_[ip].integration_weight =
69 integration_method.getWeightedPoint(ip).getWeight() *
70 sm_u.integralMeasure * sm_u.detJ;
71
72 ip_data.N_u = sm_u.N;
73 ip_data.dNdx_u = sm_u.dNdx;
74
75 // ip_data.N_p and ip_data.dNdx_p are used for both p and T variables
76 ip_data.N_p = shape_matrices[ip].N;
77 ip_data.dNdx_p = shape_matrices[ip].dNdx;
78 }
79}
80
81template <typename ShapeFunctionDisplacement, typename ShapeFunction,
82 int DisplacementDim, typename ConstitutiveTraits>
83void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
84 ShapeFunction, DisplacementDim,
85 ConstitutiveTraits>::
86 setInitialConditionsConcrete(std::vector<double> const& local_x,
87 double const t,
88 bool const /*use_monolithic_scheme*/,
89 int const /*process_id*/)
90{
91 assert(local_x.size() ==
92 temperature_size + pressure_size + displacement_size);
93
94 auto const p_L = Eigen::Map<
95 typename ShapeMatricesType::template VectorType<pressure_size> const>(
96 local_x.data() + pressure_index, pressure_size);
97
98 auto const T = Eigen::Map<typename ShapeMatricesType::template VectorType<
99 temperature_size> const>(local_x.data() + temperature_index,
100 temperature_size);
101
102 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
103 auto const& medium =
104 this->process_data_.media_map.getMedium(this->element_.getID());
105 MPL::VariableArray variables;
106
107 auto const& solid_phase = medium->phase("Solid");
108
109 unsigned const n_integration_points =
110 this->integration_method_.getNumberOfPoints();
111 for (unsigned ip = 0; ip < n_integration_points; ip++)
112 {
113 // N is used for both T and p variables.
114 auto const& N = ip_data_[ip].N_p;
115
116 ParameterLib::SpatialPosition const x_position{
117 std::nullopt, this->element_.getID(), ip,
119 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
121 this->element_, ip_data_[ip].N_u))};
122 double p_cap_ip;
123 NumLib::shapeFunctionInterpolate(-p_L, N, p_cap_ip);
124
125 variables.capillary_pressure = p_cap_ip;
126 variables.phase_pressure = -p_cap_ip;
127
128 double T_ip;
130 variables.temperature = T_ip;
131
132 this->prev_states_[ip].S_L_data->S_L =
133 medium->property(MPL::PropertyType::saturation)
134 .template value<double>(variables, x_position, t, dt);
135
136 {
137 // Set eps_m_prev from potentially non-zero eps and sigma_sw from
138 // restart.
139 SpaceTimeData const x_t{x_position, t, dt};
141 typename ConstitutiveTraits::ElasticTangentStiffnessModel{
142 this->solid_material_}
143 .eval(x_t, {T_ip, 0, {}}, C_el_data);
144
145 auto const& eps = this->current_states_[ip].eps_data.eps;
146 auto const& sigma_sw =
147 this->current_states_[ip].swelling_data.sigma_sw;
148 this->prev_states_[ip].s_mech_data->eps_m.noalias() =
149 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
150 ? eps + C_el_data.C_el.inverse() * sigma_sw
151 : eps;
152 }
153 }
154}
155
156template <typename ShapeFunctionDisplacement, typename ShapeFunction,
157 int DisplacementDim, typename ConstitutiveTraits>
158void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
159 ShapeFunction, DisplacementDim,
160 ConstitutiveTraits>::
161 assembleWithJacobian(double const t, double const dt,
162 std::vector<double> const& local_x,
163 std::vector<double> const& local_x_prev,
164 std::vector<double>& /*local_M_data*/,
165 std::vector<double>& /*local_K_data*/,
166 std::vector<double>& local_rhs_data,
167 std::vector<double>& local_Jac_data)
168{
169 auto& medium =
170 *this->process_data_.media_map.getMedium(this->element_.getID());
171
172 LocalMatrices loc_mat;
173 loc_mat.setZero();
174 LocalMatrices loc_mat_current_ip;
175 loc_mat_current_ip.setZero(); // only to set the right matrix sizes
176
177 typename ConstitutiveTraits::ConstitutiveSetting constitutive_setting;
178
179 for (unsigned ip = 0; ip < this->integration_method_.getNumberOfPoints();
180 ++ip)
181 {
182 ParameterLib::SpatialPosition const x_position{
183 std::nullopt, this->element_.getID(), ip,
185 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
187 this->element_, ip_data_[ip].N_u))};
188
189 assembleWithJacobianSingleIP(
190 t, dt, x_position, //
191 local_x, local_x_prev, //
192 ip_data_[ip], constitutive_setting,
193 medium, //
194 loc_mat_current_ip, //
195 this->current_states_[ip], this->prev_states_[ip],
196 this->material_states_[ip], this->output_data_[ip]);
197 loc_mat += loc_mat_current_ip;
198 }
199
200 massLumping(loc_mat);
201
202 addToLocalMatrixData(dt, local_x, local_x_prev, loc_mat, local_rhs_data,
203 local_Jac_data);
204}
205
206template <typename ShapeFunctionDisplacement, typename ShapeFunction,
207 int DisplacementDim, typename ConstitutiveTraits>
208void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
209 ShapeFunction, DisplacementDim,
210 ConstitutiveTraits>::
211 massLumping(typename ThermoRichardsMechanicsLocalAssembler<
212 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
213 ConstitutiveTraits>::LocalMatrices& loc_mat) const
214{
215 if (this->process_data_.apply_mass_lumping)
216 {
217 loc_mat.storage_p_a_p =
218 loc_mat.storage_p_a_p.colwise().sum().eval().asDiagonal();
219 loc_mat.storage_p_a_S =
220 loc_mat.storage_p_a_S.colwise().sum().eval().asDiagonal();
221 loc_mat.storage_p_a_S_Jpp =
222 loc_mat.storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
223 }
224}
225
226template <typename ShapeFunctionDisplacement, typename ShapeFunction,
227 int DisplacementDim, typename ConstitutiveTraits>
228void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
229 ShapeFunction, DisplacementDim,
230 ConstitutiveTraits>::
231 addToLocalMatrixData(
232 double const dt,
233 std::vector<double> const& local_x,
234 std::vector<double> const& local_x_prev,
236 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
237 ConstitutiveTraits>::LocalMatrices const& loc_mat,
238 std::vector<double>& local_rhs_data,
239 std::vector<double>& local_Jac_data) const
240{
241 constexpr auto local_matrix_dim =
242 displacement_size + pressure_size + temperature_size;
243 assert(local_x.size() == local_matrix_dim);
244
245 auto local_Jac = MathLib::createZeroedMatrix<
246 typename ShapeMatricesTypeDisplacement::template MatrixType<
247 local_matrix_dim, local_matrix_dim>>(
248 local_Jac_data, local_matrix_dim, local_matrix_dim);
249
250 auto local_rhs =
251 MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
252 template VectorType<local_matrix_dim>>(
253 local_rhs_data, local_matrix_dim);
254
255 local_Jac.noalias() = loc_mat.Jac;
256 local_rhs.noalias() = -loc_mat.res;
257
258 //
259 // -- Jacobian
260 //
261 block_TT(local_Jac).noalias() += loc_mat.M_TT / dt + loc_mat.K_TT;
262 block_Tp(local_Jac).noalias() +=
263 loc_mat.K_Tp + loc_mat.dK_TT_dp + loc_mat.M_Tp / dt;
264
265 block_pT(local_Jac).noalias() += loc_mat.M_pT / dt + loc_mat.K_pT;
266 block_pp(local_Jac).noalias() +=
267 loc_mat.K_pp + loc_mat.storage_p_a_p / dt + loc_mat.storage_p_a_S_Jpp;
268 block_pu(local_Jac).noalias() = loc_mat.M_pu / dt;
269
270 //
271 // -- Residual
272 //
273 auto const [T, p_L, u] = localDOF(local_x);
274 auto const [T_prev, p_L_prev, u_prev] = localDOF(local_x_prev);
275
276 block_T(local_rhs).noalias() -= loc_mat.M_TT * (T - T_prev) / dt +
277 loc_mat.K_TT * T + loc_mat.K_Tp * p_L +
278 loc_mat.M_Tp * (p_L - p_L_prev) / dt;
279 block_p(local_rhs).noalias() -=
280 loc_mat.K_pp * p_L + loc_mat.K_pT * T +
281 (loc_mat.storage_p_a_p + loc_mat.storage_p_a_S) * (p_L - p_L_prev) /
282 dt +
283 loc_mat.M_pu * (u - u_prev) / dt + loc_mat.M_pT * (T - T_prev) / dt;
284}
285
286template <typename ShapeFunctionDisplacement, typename ShapeFunction,
287 int DisplacementDim, typename ConstitutiveTraits>
288void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
289 ShapeFunction, DisplacementDim,
290 ConstitutiveTraits>::
291 assembleWithJacobianSingleIP(
292 double const t, double const dt,
293 ParameterLib::SpatialPosition const& x_position,
294 std::vector<double> const& local_x,
295 std::vector<double> const& local_x_prev,
297 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
298 ConstitutiveTraits>::IpData const& ip_data,
299 typename ConstitutiveTraits::ConstitutiveSetting& CS,
302 ShapeFunctionDisplacement, ShapeFunction, DisplacementDim,
303 ConstitutiveTraits>::LocalMatrices& out,
304 typename ConstitutiveTraits::StatefulData& current_state,
305 typename ConstitutiveTraits::StatefulDataPrev const& prev_state,
307 typename ConstitutiveTraits::OutputData& output_data) const
308{
309 auto const& N_u = ip_data.N_u;
310 auto const& dNdx_u = ip_data.dNdx_u;
311
312 // N and dNdx are used for both p and T variables
313 auto const& N = ip_data.N_p;
314 auto const& dNdx = ip_data.dNdx_p;
315
316 auto const x_coord =
317 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
319 this->element_, N_u);
320 auto const B =
321 LinearBMatrix::computeBMatrix<DisplacementDim,
322 ShapeFunctionDisplacement::NPOINTS,
324 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
325
326 auto const [T, p_L, u] = localDOF(local_x);
327 auto const [T_prev, p_L_prev, u_prev] = localDOF(local_x_prev);
328
329 GlobalDimVectorType const grad_T_ip = dNdx * T;
330
331 typename ConstitutiveTraits::ConstitutiveModels models(
332 this->process_data_, this->solid_material_);
333 typename ConstitutiveTraits::ConstitutiveTempData tmp;
334 typename ConstitutiveTraits::ConstitutiveData CD;
335
336 {
337 double const T_ip = N * T;
338 double const T_prev_ip = N * T_prev;
339
340 double const p_cap_ip = -N * p_L;
341 double const p_cap_prev_ip = -N * p_L_prev;
342 GlobalDimVectorType const grad_p_cap_ip = -dNdx * p_L;
343
344 KelvinVectorType eps = B * u;
345
346 CS.eval(models, t, dt, x_position, //
347 medium, //
348 {T_ip, T_prev_ip, grad_T_ip}, //
349 {p_cap_ip, p_cap_prev_ip, grad_p_cap_ip}, //
350 eps, current_state, prev_state, mat_state, tmp, output_data,
351 CD);
352 }
353
354 using NodalMatrix = typename ShapeMatricesType::NodalMatrixType;
355 NodalMatrix const NTN = N.transpose() * N;
356 NodalMatrix const dNTdN = dNdx.transpose() * dNdx;
357
358 // TODO is identity2.transpose() * B something like divergence?
359 auto const& identity2 = MathLib::KelvinVector::Invariants<
361 DisplacementDim)>::identity2;
362 typename ShapeMatricesTypeDisplacement::template MatrixType<
363 displacement_size, pressure_size>
364 BTI2N = B.transpose() * identity2 * N;
365
366 /*
367 * Conventions:
368 *
369 * * use positive signs exclusively, any negative signs must be included in
370 * the coefficients coming from the constitutive setting
371 * * the used coefficients from the constitutive setting are named after the
372 * terms they appear in, e.g. K_TT_X_dNTdN means it is a scalar (X) that
373 * will be multiplied by dNdx.transpose() * dNdx. Placefolders for the
374 * coefficients are:
375 * * X -> scalar
376 * * V -> vector
377 * * K -> Kelvin vector
378 * * the Laplace terms have a special name, e.g., K_TT_Laplace
379 * * there shall be only one contribution to each of the LocalMatrices,
380 * assigned with = assignment; this point might be relaxed in the future
381 * * this method will overwrite the data in the passed LocalMatrices& out
382 * argument, not add to it
383 */
384
385 auto const& sigma_total =
386 ConstitutiveTraits::ConstitutiveSetting::totalStress(CD, current_state);
387
388 // residual, order T, p, u
389 block_p(out.res).noalias() = dNdx.transpose() * CD.eq_p_data.rhs_p_dNT_V;
390 block_u(out.res).noalias() =
391 B.transpose() * sigma_total -
392 static_cast<int>(this->process_data_.apply_body_force_for_deformation) *
393 N_u_op(N_u).transpose() * CD.grav_data.volumetric_body_force;
394
395 // Storage matrices
396 out.storage_p_a_p.noalias() = CD.eq_p_data.storage_p_a_p_X_NTN * NTN;
397 out.storage_p_a_S.noalias() = CD.storage_data.storage_p_a_S_X_NTN * NTN;
398 out.storage_p_a_S_Jpp.noalias() =
399 CD.storage_data.storage_p_a_S_Jpp_X_NTN * NTN;
400
401 // M matrices, order T, p, u
402 out.M_TT.noalias() = CD.eq_T_data.M_TT_X_NTN * NTN;
403 out.M_Tp.noalias() = CD.vap_data.M_Tp_X_NTN * NTN;
404
405 out.M_pT.noalias() = CD.eq_p_data.M_pT_X_NTN * NTN;
406 out.M_pu.noalias() = CD.eq_p_data.M_pu_X_BTI2N * BTI2N.transpose();
407
408 // K matrices, order T, p, u
409 out.K_TT.noalias() =
410 dNdx.transpose() * CD.heat_data.K_TT_Laplace * dNdx +
411 N.transpose() * (CD.eq_T_data.K_TT_NT_V_dN.transpose() * dNdx) +
412 CD.vap_data.K_TT_X_dNTdN * dNTdN;
413
414 out.dK_TT_dp.noalias() =
415 N.transpose() * (CD.heat_data.K_Tp_NT_V_dN.transpose() * dNdx) +
416 CD.heat_data.K_Tp_X_NTN * NTN;
417 out.K_Tp.noalias() =
418 dNdx.transpose() * CD.th_osmosis_data.K_Tp_Laplace * dNdx +
419 CD.vap_data.K_Tp_X_dNTdN * dNTdN;
420
421 out.K_pp.noalias() = dNdx.transpose() * CD.eq_p_data.K_pp_Laplace * dNdx +
422 CD.vap_data.K_pp_X_dNTdN * dNTdN;
423 out.K_pT.noalias() =
424 dNdx.transpose() * CD.th_osmosis_data.K_pT_Laplace * dNdx;
425
426 // direct Jacobian contributions, order T, p, u
427 block_pT(out.Jac).noalias() = CD.vap_data.J_pT_X_dNTdN * dNTdN;
428 block_pp(out.Jac).noalias() =
429 CD.storage_data.J_pp_X_NTN * NTN +
430 CD.eq_p_data.J_pp_X_BTI2NT_u_dot_N * BTI2N.transpose() * (u - u_prev) /
431 dt * N // TODO something with volumetric strain rate?
432 + dNdx.transpose() * CD.eq_p_data.J_pp_dNT_V_N * N;
433
434 block_uT(out.Jac).noalias() =
435 B.transpose() * CD.s_mech_data.J_uT_BT_K_N * N;
436 block_up(out.Jac).noalias() =
437 B.transpose() * CD.s_mech_data.J_up_BT_K_N * N +
438 N_u_op(N_u).transpose() * CD.grav_data.J_up_HT_V_N * N;
439 block_uu(out.Jac).noalias() =
440 B.transpose() * CD.s_mech_data.stiffness_tensor * B;
441
442 out *= ip_data.integration_weight;
443}
444
445template <typename ShapeFunctionDisplacement, typename ShapeFunction,
446 int DisplacementDim, typename ConstitutiveTraits>
447void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
448 ShapeFunction, DisplacementDim,
449 ConstitutiveTraits>::
450 computeSecondaryVariableConcrete(double const t, double const dt,
451 Eigen::VectorXd const& local_x,
452 Eigen::VectorXd const& local_x_prev)
453{
454 auto const T = block_T(local_x);
455 auto const p_L = block_p(local_x);
456 auto const u = block_u(local_x);
457
458 auto const T_prev = block_T(local_x_prev);
459 auto const p_L_prev = block_p(local_x_prev);
460
461 auto const e_id = this->element_.getID();
462 auto const& process_data = this->process_data_;
463 auto& medium = *process_data.media_map.getMedium(e_id);
464
465 unsigned const n_integration_points =
466 this->integration_method_.getNumberOfPoints();
467
468 double saturation_avg = 0;
469 double porosity_avg = 0;
470 double liquid_density_avg = 0;
471 double viscosity_avg = 0;
472
474 KV sigma_avg = KV::Zero();
475
476 typename ConstitutiveTraits::ConstitutiveSetting constitutive_setting;
477
478 typename ConstitutiveTraits::ConstitutiveModels models(
479 process_data, this->solid_material_);
480 typename ConstitutiveTraits::ConstitutiveTempData tmp;
481 typename ConstitutiveTraits::ConstitutiveData CD;
482
483 for (unsigned ip = 0; ip < n_integration_points; ip++)
484 {
485 auto& current_state = this->current_states_[ip];
486 auto& output_data = this->output_data_[ip];
487
488 auto const& ip_data = ip_data_[ip];
489
490 // N is used for both p and T variables
491 auto const& N = ip_data.N_p;
492 auto const& N_u = ip_data.N_u;
493 auto const& dNdx_u = ip_data.dNdx_u;
494 auto const& dNdx = ip_data.dNdx_p;
495
496 ParameterLib::SpatialPosition const x_position{
497 std::nullopt, this->element_.getID(), ip,
499 NumLib::interpolateCoordinates<ShapeFunctionDisplacement,
501 this->element_, N_u))};
502 auto const x_coord =
503 NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
505 this->element_, N_u);
506 auto const B =
507 LinearBMatrix::computeBMatrix<DisplacementDim,
508 ShapeFunctionDisplacement::NPOINTS,
510 dNdx_u, N_u, x_coord, this->is_axially_symmetric_);
511
512 double const T_ip = N * T;
513 double const T_prev_ip = N * T_prev;
514 GlobalDimVectorType const grad_T_ip = dNdx * T;
515
516 double const p_cap_ip = -N * p_L;
517 double const p_cap_prev_ip = -N * p_L_prev;
518 GlobalDimVectorType const grad_p_cap_ip = -dNdx * p_L;
519
520 KelvinVectorType eps = B * u;
521
522 constitutive_setting.eval(models, //
523 t, dt, x_position, //
524 medium, //
525 {T_ip, T_prev_ip, grad_T_ip}, //
526 {p_cap_ip, p_cap_prev_ip, grad_p_cap_ip}, //
527 eps, current_state, this->prev_states_[ip],
528 this->material_states_[ip], tmp, output_data,
529 CD);
530
531 saturation_avg += current_state.S_L_data.S_L;
532 porosity_avg += current_state.poro_data.phi;
533
534 liquid_density_avg += output_data.rho_L_data.rho_LR;
535 viscosity_avg += output_data.mu_L_data.viscosity;
536 sigma_avg += ConstitutiveTraits::ConstitutiveSetting::statefulStress(
537 current_state);
538 }
539 saturation_avg /= n_integration_points;
540 porosity_avg /= n_integration_points;
541 viscosity_avg /= n_integration_points;
542 liquid_density_avg /= n_integration_points;
543 sigma_avg /= n_integration_points;
544
545 (*process_data.element_saturation)[e_id] = saturation_avg;
546 (*process_data.element_porosity)[e_id] = porosity_avg;
547 (*process_data.element_liquid_density)[e_id] = liquid_density_avg;
548 (*process_data.element_viscosity)[e_id] = viscosity_avg;
549
550 Eigen::Map<KV>(
551 &(*process_data.element_stresses)[e_id * KV::RowsAtCompileTime]) =
553
555 ShapeFunction, typename ShapeFunctionDisplacement::MeshElement,
556 DisplacementDim>(this->element_, this->is_axially_symmetric_, p_L,
557 *process_data.pressure_interpolated);
559 ShapeFunction, typename ShapeFunctionDisplacement::MeshElement,
560 DisplacementDim>(this->element_, this->is_axially_symmetric_, T,
561 *process_data.temperature_interpolated);
562}
563} // namespace ThermoRichardsMechanics
564} // namespace ProcessLib
double getWeight() const
Definition: WeightedPoint.h:80
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
MatrixType< _kelvin_vector_size, _number_of_dof > BMatrixType
Definition: BMatrixPolicy.h:55
ShapeMatrixPolicyType< ShapeFunctionDisplacement, DisplacementDim > ShapeMatricesTypeDisplacement
ShapeMatrixPolicyType< ShapeFunction, DisplacementDim > ShapeMatricesType
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
Eigen::Matrix< double, 4, 1 > kelvinVectorToSymmetricTensor(Eigen::Matrix< double, 4, 1, Eigen::ColMajor, 4, 1 > const &v)
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:24
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)
Definition: EigenMapTools.h:32
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition: Interpolation.h:27
double interpolateXCoordinate(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
void interpolateToHigherOrderNodes(MeshLib::Element const &element, bool const is_axially_symmetric, Eigen::MatrixBase< EigenMatrixType > const &node_values, MeshLib::PropertyVector< double > &interpolated_values_global_vector)
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)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
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.
Definition: LinearBMatrix.h:42
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType