103 ShapeFunctionPressure, DisplacementDim>::
104 assembleWithJacobian(
double const t,
double const dt,
105 std::vector<double>
const& local_x,
106 std::vector<double>
const& local_x_prev,
107 std::vector<double>& local_rhs_data,
108 std::vector<double>& local_Jac_data)
110 assert(local_x.size() == pressure_size + displacement_size);
112 auto p = Eigen::Map<
typename ShapeMatricesTypePressure::template VectorType<
113 pressure_size>
const>(local_x.data() + pressure_index, pressure_size);
116 Eigen::Map<
typename ShapeMatricesTypeDisplacement::template VectorType<
117 displacement_size>
const>(local_x.data() + displacement_index,
121 Eigen::Map<
typename ShapeMatricesTypePressure::template VectorType<
122 pressure_size>
const>(local_x_prev.data() + pressure_index,
125 Eigen::Map<
typename ShapeMatricesTypeDisplacement::template VectorType<
126 displacement_size>
const>(local_x_prev.data() + displacement_index,
130 typename ShapeMatricesTypeDisplacement::template MatrixType<
131 displacement_size + pressure_size,
132 displacement_size + pressure_size>>(
133 local_Jac_data, displacement_size + pressure_size,
134 displacement_size + pressure_size);
137 typename ShapeMatricesTypeDisplacement::template VectorType<
138 displacement_size + pressure_size>>(
139 local_rhs_data, displacement_size + pressure_size);
142 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
146 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
150 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
153 typename ShapeMatricesTypeDisplacement::template MatrixType<
154 displacement_size, pressure_size>
155 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
156 displacement_size, pressure_size>::Zero(displacement_size,
159 typename ShapeMatricesTypeDisplacement::template MatrixType<
160 pressure_size, displacement_size>
161 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
162 pressure_size, displacement_size>::Zero(pressure_size,
165 typename ShapeMatricesTypeDisplacement::template MatrixType<
166 pressure_size, displacement_size>
167 Kpu_k = ShapeMatricesTypeDisplacement::template MatrixType<
168 pressure_size, displacement_size>::Zero(pressure_size,
171 auto const& solid_material =
173 _process_data.solid_materials, _process_data.material_ids,
179 unsigned const n_integration_points =
180 _integration_method.getNumberOfPoints();
182 auto const& b = _process_data.specific_body_force;
183 auto const& medium = _process_data.media_map.getMedium(_element.getID());
184 auto const& solid = medium->phase(
"Solid");
185 auto const& fluid = fluidPhase(*medium);
191 medium->property(MPL::PropertyType::reference_temperature)
192 .template value<double>(vars, x_position, t, dt);
195 auto const& identity2 = Invariants::identity2;
197 for (
unsigned ip = 0; ip < n_integration_points; ip++)
199 auto const& w = _ip_data[ip].integration_weight;
201 auto const& N_u = _ip_data[ip].N_u;
202 auto const& dNdx_u = _ip_data[ip].dNdx_u;
204 auto const& N_p = _ip_data[ip].N_p;
205 auto const& dNdx_p = _ip_data[ip].dNdx_p;
208 std::nullopt, _element.getID(),
218 ShapeFunctionDisplacement::NPOINTS,
220 dNdx_u, N_u, x_coord, _is_axially_symmetric);
222 auto& eps = _ip_data[ip].eps;
223 eps.noalias() = B * u;
224 auto const& sigma_eff = _ip_data[ip].sigma_eff;
226 double const p_int_pt = N_p.dot(p);
228 phase_pressure = p_int_pt;
230 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
231 t, x_position, dt, T_ref);
232 auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
234 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
235 .template value<double>(vars, x_position, t, dt);
238 solid.property(MPL::PropertyType::density)
239 .template value<double>(vars, x_position, t, dt);
240 auto const porosity =
241 medium->property(MPL::PropertyType::porosity)
242 .template value<double>(vars, x_position, t, dt);
244 auto const [rho_fr, mu] =
249 fluid.property(MPL::PropertyType::density)
250 .template dValue<double>(vars, _process_data.phase_variable,
257 auto const sigma_total =
258 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
267 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
273 medium->property(MPL::PropertyType::permeability)
274 .value(vars, x_position, t, dt));
277 (*medium)[MPL::PropertyType::permeability].dValue(
278 vars, MPL::Variable::mechanical_strain, x_position, t, dt));
280 auto const K_over_mu = K / mu;
282 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
289 .template block<displacement_size, displacement_size>(
290 displacement_index, displacement_index)
291 .noalias() += B.transpose() * C * B * w;
293 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
294 local_rhs.template segment<displacement_size>(displacement_index)
296 (B.transpose() * sigma_eff - N_u_op(N_u).transpose() * rho * b) * w;
301 Kup.noalias() += B.transpose() * alpha * identity2 * N_p * w;
306 laplace_p.noalias() +=
307 rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
309 storage_p.noalias() +=
310 rho_fr * N_p.transpose() * N_p * w *
311 ((alpha - porosity) * (1.0 - alpha) / K_S + porosity * beta_p);
315 add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
317 (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
319 local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
320 dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
326 rho_fr * alpha * N_p.transpose() * identity2.transpose() * B * w;
331 dNdx_p * p - rho_fr * b) *
332 dkde * B * rho_fr / mu * w;
336 .template block<displacement_size, pressure_size>(displacement_index,
340 if (_process_data.mass_lumping)
342 storage_p = storage_p.colwise().sum().eval().asDiagonal();
344 if constexpr (pressure_size == displacement_size)
346 Kpu = Kpu.colwise().sum().eval().asDiagonal();
347 Kpu_k = Kpu_k.colwise().sum().eval().asDiagonal();
353 .template block<pressure_size, pressure_size>(pressure_index,
355 .noalias() += laplace_p + storage_p / dt + add_p_derivative;
359 .template block<pressure_size, displacement_size>(pressure_index,
361 .noalias() += Kpu / dt + Kpu_k;
364 local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
365 laplace_p * p + storage_p * (p - p_prev) / dt + Kpu * (u - u_prev) / dt;
368 local_rhs.template segment<displacement_size>(displacement_index)
369 .noalias() += Kup * p;
375 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
376 getIntPtDarcyVelocity(
378 std::vector<GlobalVector*>
const& x,
379 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
380 std::vector<double>& cache)
const
382 int const hydraulic_process_id = _process_data.hydraulic_process_id;
385 assert(!indices.empty());
386 auto const local_x = x[hydraulic_process_id]->get(indices);
388 unsigned const n_integration_points =
389 _integration_method.getNumberOfPoints();
392 double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
393 cache, DisplacementDim, n_integration_points);
395 auto p = Eigen::Map<
typename ShapeMatricesTypePressure::template VectorType<
396 pressure_size>
const>(local_x.data() + pressure_index, pressure_size);
401 auto const& medium = _process_data.media_map.getMedium(_element.getID());
402 auto const& fluid = fluidPhase(*medium);
409 double const dt = std::numeric_limits<double>::quiet_NaN();
411 medium->property(MPL::PropertyType::reference_temperature)
412 .template value<double>(vars, x_position, t, dt);
414 auto const& identity2 = Invariants::identity2;
416 for (
unsigned ip = 0; ip < n_integration_points; ip++)
419 std::nullopt, _element.getID(),
423 _element, _ip_data[ip].N_u))};
425 double const p_int_pt = _ip_data[ip].N_p.dot(p);
427 phase_pressure = p_int_pt;
429 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
430 .template value<double>(vars, x_position, t, dt);
434 auto const sigma_total =
435 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
442 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
448 medium->property(MPL::PropertyType::permeability)
449 .value(vars, x_position, t, dt));
451 auto const [rho_fr, mu] =
455 auto const K_over_mu = K / mu;
457 auto const& b = _process_data.specific_body_force;
460 auto const& dNdx_p = _ip_data[ip].dNdx_p;
461 cache_matrix.col(ip).noalias() =
462 -K_over_mu * dNdx_p * p + K_over_mu * rho_fr * b;
471 ShapeFunctionPressure, DisplacementDim>::
472 assembleWithJacobianForPressureEquations(
473 const double t,
double const dt, Eigen::VectorXd
const& local_x,
474 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
475 std::vector<double>& local_Jac_data)
479 template VectorType<pressure_size>>(
480 local_b_data, pressure_size);
485 auto const p = local_x.template segment<pressure_size>(pressure_index);
488 local_x_prev.template segment<pressure_size>(pressure_index);
491 typename ShapeMatricesTypeDisplacement::template MatrixType<
492 pressure_size, pressure_size>>(local_Jac_data, pressure_size,
496 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
500 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
504 ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
507 auto const& solid_material =
509 _process_data.solid_materials, _process_data.material_ids,
515 auto const& medium = _process_data.media_map.getMedium(_element.getID());
516 auto const& fluid = fluidPhase(*medium);
522 medium->property(MPL::PropertyType::reference_temperature)
523 .template value<double>(vars, x_position, t, dt);
526 auto const& identity2 = Invariants::identity2;
528 auto const staggered_scheme =
529 std::get<Staggered>(_process_data.coupling_scheme);
530 auto const fixed_stress_stabilization_parameter =
531 staggered_scheme.fixed_stress_stabilization_parameter;
532 auto const fixed_stress_over_time_step =
533 staggered_scheme.fixed_stress_over_time_step;
535 int const n_integration_points = _integration_method.getNumberOfPoints();
536 for (
int ip = 0; ip < n_integration_points; ip++)
538 auto const& w = _ip_data[ip].integration_weight;
540 auto const& N_p = _ip_data[ip].N_p;
541 auto const& dNdx_p = _ip_data[ip].dNdx_p;
544 std::nullopt, _element.getID(),
548 _element, _ip_data[ip].N_p))};
550 double const p_int_pt = N_p.dot(p);
552 phase_pressure = p_int_pt;
554 auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
555 t, x_position, dt, T_ref);
556 auto const K_S = solid_material.getBulkModulus(t, x_position, &C_el);
559 medium->property(MPL::PropertyType::biot_coefficient)
560 .template value<double>(vars, x_position, t, dt);
564 auto const sigma_total =
565 (_ip_data[ip].sigma_eff - alpha_b * p_int_pt * identity2).eval();
572 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
578 medium->property(MPL::PropertyType::permeability)
579 .value(vars, x_position, t, dt));
580 auto const porosity =
581 medium->property(MPL::PropertyType::porosity)
582 .template value<double>(vars, x_position, t, dt);
584 auto const [rho_fr, mu] =
589 fluid.property(MPL::PropertyType::density)
590 .template dValue<double>(vars, _process_data.phase_variable,
594 auto const K_over_mu = K / mu;
597 rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
601 fixed_stress_stabilization_parameter * alpha_b * alpha_b / K_S;
603 storage.noalias() += rho_fr * N_p.transpose() * N_p * w *
604 ((alpha_b - porosity) * (1.0 - alpha_b) / K_S +
605 porosity * beta_p + beta_FS);
607 auto const& b = _process_data.specific_body_force;
610 local_rhs.noalias() +=
611 dNdx_p.transpose() * rho_fr * rho_fr * K_over_mu * b * w;
616 add_p_derivative.noalias() += rho_fr * beta_p * dNdx_p.transpose() *
618 (dNdx_p * p - 2.0 * rho_fr * b) * N_p * w;
620 if (!fixed_stress_over_time_step)
622 auto const& eps = _ip_data[ip].eps;
623 auto const& eps_prev = _ip_data[ip].eps_prev;
624 const double eps_v_dot =
625 (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
628 double const strain_rate_b =
629 alpha_b * eps_v_dot -
630 beta_FS * _ip_data[ip].strain_rate_variable;
632 local_rhs.noalias() -= strain_rate_b * rho_fr * N_p * w;
637 local_rhs.noalias() -=
638 alpha_b * _ip_data[ip].strain_rate_variable * rho_fr * N_p * w;
641 local_Jac.noalias() = laplace + storage / dt + add_p_derivative;
643 local_rhs.noalias() -= laplace * p + storage * (p - p_prev) / dt;
649 ShapeFunctionPressure, DisplacementDim>::
650 assembleWithJacobianForDeformationEquations(
651 const double t,
double const dt, Eigen::VectorXd
const& local_x,
652 std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
654 auto const p = local_x.template segment<pressure_size>(pressure_index);
656 local_x.template segment<displacement_size>(displacement_index);
659 typename ShapeMatricesTypeDisplacement::template MatrixType<
660 displacement_size, displacement_size>>(
661 local_Jac_data, displacement_size, displacement_size);
665 template VectorType<displacement_size>>(
666 local_b_data, displacement_size);
671 auto const& medium = _process_data.media_map.getMedium(_element.getID());
672 auto const& solid = medium->phase(
"Solid");
673 auto const& fluid = fluidPhase(*medium);
679 medium->property(MPL::PropertyType::reference_temperature)
680 .template value<double>(vars, x_position, t, dt);
683 int const n_integration_points = _integration_method.getNumberOfPoints();
684 for (
int ip = 0; ip < n_integration_points; ip++)
686 auto const& w = _ip_data[ip].integration_weight;
688 auto const& N_u = _ip_data[ip].N_u;
689 auto const& dNdx_u = _ip_data[ip].dNdx_u;
691 auto const& N_p = _ip_data[ip].N_p;
694 std::nullopt, _element.getID(),
700 auto const x_coord = x_position.getCoordinates().value()[0];
703 ShapeFunctionDisplacement::NPOINTS,
705 dNdx_u, N_u, x_coord, _is_axially_symmetric);
707 auto& eps = _ip_data[ip].eps;
708 auto const& sigma_eff = _ip_data[ip].sigma_eff;
710 phase_pressure = N_p.dot(p);
712 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
713 .template value<double>(vars, x_position, t, dt);
715 solid.property(MPL::PropertyType::density)
716 .template value<double>(vars, x_position, t, dt);
717 auto const porosity =
718 medium->property(MPL::PropertyType::porosity)
719 .template value<double>(vars, x_position, t, dt);
722 t, dt, x_position, fluid, vars);
724 auto const& b = _process_data.specific_body_force;
727 DisplacementDim)>::identity2;
729 eps.noalias() = B * u;
734 auto C = _ip_data[ip].updateConstitutiveRelation(vars, t, x_position,
737 local_Jac.noalias() += B.transpose() * C * B * w;
742 double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
743 local_rhs.noalias() -=
744 (B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) -
745 N_u_op(N_u).transpose() * rho * b) *
777 ShapeFunctionPressure, DisplacementDim>::
778 setInitialConditionsConcrete(Eigen::VectorXd
const local_x,
785 auto const& medium = _process_data.media_map.getMedium(_element.getID());
787 auto const p = local_x.template segment<pressure_size>(pressure_index);
789 local_x.template segment<displacement_size>(displacement_index);
791 auto const& identity2 = Invariants::identity2;
792 const double dt = 0.0;
796 int const n_integration_points = _integration_method.getNumberOfPoints();
797 for (
int ip = 0; ip < n_integration_points; ip++)
799 auto const& N_u = _ip_data[ip].N_u;
800 auto const& dNdx_u = _ip_data[ip].dNdx_u;
803 std::nullopt, _element.getID(),
812 ShapeFunctionDisplacement::NPOINTS,
814 dNdx_u, N_u, x_coord, _is_axially_symmetric);
816 auto& eps = _ip_data[ip].eps;
817 eps.noalias() = B * u;
822 if (_process_data.initial_stress.isTotalStress())
824 auto const& N_p = _ip_data[ip].N_p;
826 medium->property(MPL::PropertyType::biot_coefficient)
827 .template value<double>(vars, x_position, t, dt);
829 auto& sigma_eff = _ip_data[ip].sigma_eff;
830 sigma_eff.noalias() += alpha_b * N_p.dot(p) * identity2;
831 _ip_data[ip].sigma_eff_prev.noalias() = sigma_eff;
839 ShapeFunctionPressure, DisplacementDim>::
840 postNonLinearSolverConcrete(Eigen::VectorXd
const& local_x,
841 Eigen::VectorXd
const& local_x_prev,
842 double const t,
double const dt,
843 int const process_id)
849 int const n_integration_points = _integration_method.getNumberOfPoints();
851 auto const staggered_scheme_ptr =
852 std::get_if<Staggered>(&_process_data.coupling_scheme);
854 if (staggered_scheme_ptr &&
855 process_id == _process_data.hydraulic_process_id)
857 if (!staggered_scheme_ptr->fixed_stress_over_time_step)
860 local_x.template segment<pressure_size>(pressure_index);
862 local_x_prev.template segment<pressure_size>(pressure_index);
864 for (
int ip = 0; ip < n_integration_points; ip++)
866 auto& ip_data = _ip_data[ip];
868 auto const& N_p = ip_data.N_p;
870 ip_data.strain_rate_variable = N_p.dot(p - p_prev) / dt;
875 if (!staggered_scheme_ptr ||
876 process_id == _process_data.mechanics_related_process_id)
882 _process_data.media_map.getMedium(_element.getID());
885 medium->property(MPL::PropertyType::reference_temperature)
890 local_x.template segment<displacement_size>(displacement_index);
895 for (
int ip = 0; ip < n_integration_points; ip++)
897 auto const& N_u = _ip_data[ip].N_u;
898 auto const& dNdx_u = _ip_data[ip].dNdx_u;
900 x_position = {std::nullopt, _element.getID(),
902 ShapeFunctionDisplacement,
908 DisplacementDim, ShapeFunctionDisplacement::NPOINTS,
910 _is_axially_symmetric);
912 auto& eps = _ip_data[ip].eps;
913 eps.noalias() = B * u;
917 _ip_data[ip].updateConstitutiveRelation(vars, t, x_position, dt, u,
927 DisplacementDim>::postTimestepConcrete(Eigen::VectorXd
const& local_x,
928 Eigen::VectorXd
const& local_x_prev,
929 double const t,
double const dt,
930 int const process_id)
932 auto const staggered_scheme_ptr =
933 std::get_if<Staggered>(&_process_data.coupling_scheme);
935 if (staggered_scheme_ptr &&
936 process_id == _process_data.hydraulic_process_id)
938 if (staggered_scheme_ptr->fixed_stress_over_time_step)
940 auto const fixed_stress_stabilization_parameter =
941 staggered_scheme_ptr->fixed_stress_stabilization_parameter;
944 local_x.template segment<pressure_size>(pressure_index);
946 local_x_prev.template segment<pressure_size>(pressure_index);
951 auto const& solid_material =
953 _process_data.solid_materials, _process_data.material_ids,
957 _process_data.media_map.getMedium(_element.getID());
961 medium->property(MPL::PropertyType::reference_temperature)
962 .template value<double>(vars, x_position, t, dt);
965 int const n_integration_points =
966 _integration_method.getNumberOfPoints();
967 for (
int ip = 0; ip < n_integration_points; ip++)
969 auto& ip_data = _ip_data[ip];
971 auto const& N_p = ip_data.N_p;
974 std::nullopt, _element.getID(),
980 auto const& eps = ip_data.eps;
981 auto const& eps_prev = ip_data.eps_prev;
982 const double eps_v_dot =
983 (Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
985 auto const C_el = ip_data.computeElasticTangentStiffness(
986 t, x_position, dt, T_ref);
988 solid_material.getBulkModulus(t, x_position, &C_el);
991 medium->property(MPL::PropertyType::biot_coefficient)
992 .template value<double>(vars, x_position, t, dt);
994 ip_data.strain_rate_variable =
995 eps_v_dot - fixed_stress_stabilization_parameter * alpha_b *
996 N_p.dot(p - p_prev) / dt / K_S;
1001 unsigned const n_integration_points =
1002 _integration_method.getNumberOfPoints();
1004 for (
unsigned ip = 0; ip < n_integration_points; ip++)
1006 _ip_data[ip].pushBackState();
1116 ShapeFunctionPressure, DisplacementDim>::
1117 computeSecondaryVariableConcrete(
double const t,
double const dt,
1118 Eigen::VectorXd
const& local_x,
1119 Eigen::VectorXd
const& )
1121 auto const p = local_x.template segment<pressure_size>(pressure_index);
1124 ShapeFunctionPressure,
typename ShapeFunctionDisplacement::MeshElement,
1125 DisplacementDim>(_element, _is_axially_symmetric, p,
1126 *_process_data.pressure_interpolated);
1128 int const elem_id = _element.getID();
1130 unsigned const n_integration_points =
1131 _integration_method.getNumberOfPoints();
1133 auto const& medium = _process_data.media_map.getMedium(elem_id);
1140 Eigen::Matrix<double, 3, 3>::Zero());
1142 auto const& identity2 = Invariants::identity2;
1144 for (
unsigned ip = 0; ip < n_integration_points; ip++)
1146 auto const& eps = _ip_data[ip].eps;
1147 sigma_eff_sum += _ip_data[ip].sigma_eff;
1150 std::nullopt, _element.getID(),
1154 _element, _ip_data[ip].N_u))};
1156 auto const alpha = medium->property(MPL::PropertyType::biot_coefficient)
1157 .template value<double>(vars, x_position, t, dt);
1158 double const p_int_pt = _ip_data[ip].N_p.dot(p);
1160 phase_pressure = p_int_pt;
1165 auto const sigma_total =
1166 (_ip_data[ip].sigma_eff - alpha * p_int_pt * identity2).eval();
1174 _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
1180 medium->property(MPL::PropertyType::permeability)
1181 .value(vars, x_position, t, dt));
1184 Eigen::Map<Eigen::VectorXd>(
1185 &(*_process_data.permeability)[elem_id * KelvinVectorSize],
1186 KelvinVectorSize) = k_sum / n_integration_points;
1188 Eigen::Matrix<double, 3, 3, 0, 3, 3>
const sigma_avg =
1190 n_integration_points;
1192 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(sigma_avg);
1194 Eigen::Map<Eigen::Vector3d>(
1195 &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
1198 auto eigen_vectors = e_s.eigenvectors();
1200 for (
auto i = 0; i < 3; i++)
1202 Eigen::Map<Eigen::Vector3d>(
1203 &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
1204 eigen_vectors.col(i);