42 std::span<double const, 3>
const nat_coords)
44 using ShapeMatricesType =
47 typename ShapeMatricesType::ShapeMatrices shape_matrix(
56 using ElementTrait = std::remove_pointer_t<
decltype(tag)>;
57 using ShapeFunction =
typename ElementTrait::ShapeFunction;
58 using MeshElement =
typename ElementTrait::Element;
61 dynamic_cast<MeshElement const*
>(&bulk_element) !=
nullptr)
64 ShapeFunction, ShapeMatricesType>(bulk_element);
66 fe.template computeShapeFunctions<NumLib::ShapeMatrixType::N>(
67 nat_coords.data(), shape_matrix, GlobalDim,
76 OGS_FATAL(
"Element type '{:s}' is not supported as anchor element.",
80 return shape_matrix.N;
85 Eigen::Vector<double, GlobalDim>
const& f,
86 Eigen::Matrix<double, GlobalDim, GlobalDim>
const& Df,
87 std::vector<Eigen::RowVectorXd>
const& shape_matrices,
88 std::size_t
const num_dof,
89 std::array<std::size_t, 2>
const& nodes_per_element)
92 constexpr auto even_odd_sign = [](std::size_t
const n)
93 {
return (n % 2 == 0) ? 1.0 : -1.0; };
95 Eigen::VectorXd local_rhs(num_dof);
96 Eigen::MatrixXd local_Jac(num_dof, num_dof);
100 for (std::size_t node_idx = 0; node_idx < 2; ++node_idx)
105 local_rhs.segment(node_idx * GlobalDim * nodes_per_element[0],
106 GlobalDim * nodes_per_element[node_idx]) =
107 even_odd_sign(node_idx) *
108 shape_matrices[node_idx]
110 .replicate<GlobalDim, 1>()
111 .cwiseProduct(f.transpose()
112 .replicate(nodes_per_element[node_idx], 1)
114 for (std::size_t node_idx_inner = 0; node_idx_inner < 2;
117 Eigen::MatrixXd
const& ones = Eigen::MatrixXd::Ones(
118 nodes_per_element[node_idx], nodes_per_element[node_idx_inner]);
124 local_Jac.block(node_idx * GlobalDim * nodes_per_element[0],
125 node_idx_inner * GlobalDim * nodes_per_element[0],
126 GlobalDim * nodes_per_element[node_idx],
127 GlobalDim * nodes_per_element[node_idx_inner]) =
128 (even_odd_sign(node_idx) *
129 shape_matrices[node_idx]
131 .replicate<GlobalDim, 1>() *
132 even_odd_sign(node_idx_inner) *
133 shape_matrices[node_idx_inner].replicate<1, GlobalDim>())
134 .cwiseProduct(kroneckerProduct(Df, ones));
137 return {local_rhs, local_Jac};
144 std::size_t
const source_term_mesh_id,
146 const int variable_id)
147 : dof_table_bulk_(dof_table_bulk),
148 bulk_mesh_(bulk_mesh),
149 source_term_mesh_id_(source_term_mesh_id),
151 variable_id_(variable_id),
155 std::array<int, GlobalDim> arr{};
156 std::iota(arr.begin(), arr.end(), 0);
160 DBUG(
"Create EmbeddedAnchor.");
161 std::string_view
const bulk_element_ids_string =
"bulk_element_ids";
163 st_mesh_.getProperties().template getPropertyVector<std::size_t>(
164 bulk_element_ids_string);
166 std::string_view
const natural_coordinates_string =
"natural_coordinates";
167 natural_coordinates_ =
168 st_mesh_.getProperties().template getPropertyVector<double>(
169 natural_coordinates_string);
171 std::string_view
const maximum_anchor_stress_string =
172 "maximum_anchor_stress";
173 maximum_anchor_stress_ =
174 st_mesh_.getProperties().template getPropertyVector<double>(
175 maximum_anchor_stress_string);
177 std::string_view
const initial_anchor_stress_string =
178 "initial_anchor_stress";
179 initial_anchor_stress_ =
180 st_mesh_.getProperties().template getPropertyVector<double>(
181 initial_anchor_stress_string);
183 std::string_view
const residual_anchor_stress_string =
184 "residual_anchor_stress";
185 residual_anchor_stress_ =
186 st_mesh_.getProperties().template getPropertyVector<double>(
187 residual_anchor_stress_string);
189 std::string_view
const anchor_cross_sectional_area_string =
190 "anchor_cross_sectional_area";
191 cross_sectional_area_ =
192 st_mesh_.getProperties().template getPropertyVector<double>(
193 anchor_cross_sectional_area_string);
195 std::string_view
const anchor_stiffness_string =
"anchor_stiffness";
197 st_mesh_.getProperties().template getPropertyVector<double>(
198 anchor_stiffness_string);
205 std::array<std::size_t, 2>& nodes_per_element,
206 std::vector<Eigen::RowVectorXd>& shape_matrices,
207 std::vector<GlobalIndexType>& global_indices,
208 Eigen::Vector<double, 2 * GlobalDim>& local_x,
212 for (
auto&& [node_idx, anchor_node_id] : ranges::views::enumerate(
217 auto const bulk_element_id = (*bulk_element_ids_)[anchor_node_id];
220 auto const& bulk_element = *bulk_mesh_.getElement(bulk_element_id);
221 nodes_per_element[node_idx] = bulk_element.nodes().size();
223 std::span<double const, 3>
const nat_coords(
224 &(*natural_coordinates_)[3 * anchor_node_id], 3);
226 shape_matrices.push_back(std::move(N));
228 for (
int component = 0; component < GlobalDim; ++component)
230 auto const& global_indices_percomponent = getIndices(
231 variable_id_, component, bulk_element_id, dof_table_bulk_);
233 Eigen::VectorXd
const local_x_element_percomponent =
236 global_indices.insert(global_indices.cend(),
237 global_indices_percomponent.begin(),
238 global_indices_percomponent.end());
245 local_x_element_percomponent,
246 shape_matrices[node_idx],
258 DBUG(
"Assemble EmbeddedAnchor.");
260 using GlobalDimVector = Eigen::Vector<double, GlobalDim>;
261 using GlobalDimMatrix = Eigen::Matrix<double, GlobalDim, GlobalDim>;
265 auto const anchor_element_id = anchor_element->getID();
266 std::vector<GlobalIndexType> global_indices;
267 Eigen::Vector<double, 2 * GlobalDim> local_x;
268 std::vector<Eigen::RowVectorXd> shape_matrices;
269 std::array<std::size_t, 2> nodes_per_element;
271 getShapeMatricesAndGlobalIndicesAndDisplacements(
272 anchor_element, nodes_per_element, shape_matrices, global_indices,
275 auto node_coords = [anchor_element](
int const i)
276 {
return anchor_element->getNode(i)->asEigenVector3d(); };
277 GlobalDimVector
const l_original =
278 (node_coords(1) - node_coords(0)).
template head<GlobalDim>();
279 double const l_original_norm = l_original.norm();
282 auto u = [&local_x](
int const i)
284 GlobalDimVector
const l = l_original + u(1) - u(0);
286 double const K = (*cross_sectional_area_)[anchor_element_id] *
287 (*anchor_stiffness_)[anchor_element_id];
288 double const initial_force =
289 (*cross_sectional_area_)[anchor_element_id] *
290 (*initial_anchor_stress_)[anchor_element_id];
291 double const max_force = (*cross_sectional_area_)[anchor_element_id] *
292 (*maximum_anchor_stress_)[anchor_element_id];
293 double const residual_force =
294 (*cross_sectional_area_)[anchor_element_id] *
295 (*residual_anchor_stress_)[anchor_element_id];
297 double const strain = (l.norm() - l_original_norm) / l_original_norm;
299 GlobalDimVector
const f_friction =
300 residual_force * l_original / l_original_norm;
301 GlobalDimVector
const f_elastic =
302 l_original / l_original_norm * (initial_force + K * strain);
303 GlobalDimVector
const f =
304 ((f_elastic.norm() < max_force)) ? f_elastic : f_friction;
306 GlobalDimMatrix
const Df_friction = GlobalDimMatrix::Zero();
307 GlobalDimMatrix
const Df_elastic = l_original / l_original_norm * K *
308 l.transpose() / l.norm() /
310 GlobalDimMatrix
const Df =
311 (f_elastic.norm() < max_force) ? Df_elastic : Df_friction;
314 f, Df, shape_matrices, global_indices.size(), nodes_per_element);
316 b.
add(global_indices, local_rhs);
319 jac->
add({global_indices, global_indices}, local_Jac);
void getShapeMatricesAndGlobalIndicesAndDisplacements(MeshLib::Element const *const anchor_element, std::array< std::size_t, 2 > &nodes_per_element, std::vector< Eigen::RowVectorXd > &shape_matrices, std::vector< GlobalIndexType > &global_indices, Eigen::Vector< double, 2 *GlobalDim > &local_x, GlobalVector const &x, ParameterLib::SpatialPosition &pos) const
std::tuple< Eigen::VectorXd, Eigen::MatrixXd > assembleLocalBJac(Eigen::Vector< double, GlobalDim > const &f, Eigen::Matrix< double, GlobalDim, GlobalDim > const &Df, std::vector< Eigen::RowVectorXd > const &shape_matrices, std::size_t const num_dof, std::array< std::size_t, 2 > const &nodes_per_element)