OGS
EmbeddedAnchor.cpp
Go to the documentation of this file.
1
11#include "EmbeddedAnchor.h"
12
13#include <Eigen/Core>
14#include <cassert>
15#include <numeric>
16#include <range/v3/view/enumerate.hpp>
17#include <unsupported/Eigen/KroneckerProduct>
18
28
29namespace ProcessLib
30{
31// The local indices are, due to the nature of the DOF table, all even
32// for the first node and odd for the second node.
33template <int GlobalDim>
34auto nodeLocalIndices(std::size_t const i)
35{
36 return Eigen::seqN(i, Eigen::fix<GlobalDim>, Eigen::fix<2>);
37}
38
39template <int GlobalDim>
40Eigen::RowVectorXd computeShapeMatrix(
41 MeshLib::Element const& bulk_element,
42 std::span<double const, 3> const nat_coords)
43{
44 using ShapeMatricesType =
45 EigenDynamicShapeMatrixPolicy<void /* dummy */, -1 /* dynamic order */>;
46
47 typename ShapeMatricesType::ShapeMatrices shape_matrix(
48 bulk_element.getDimension(), GlobalDim,
49 bulk_element.getNumberOfNodes());
50
51 bool matched = false;
52
54 [&](auto* tag)
55 {
56 using ElementTrait = std::remove_pointer_t<decltype(tag)>;
57 using ShapeFunction = typename ElementTrait::ShapeFunction;
58 using MeshElement = typename ElementTrait::Element;
59
60 if (!matched &&
61 dynamic_cast<MeshElement const*>(&bulk_element) != nullptr)
62 {
64 ShapeFunction, ShapeMatricesType>(bulk_element);
65
66 fe.template computeShapeFunctions<NumLib::ShapeMatrixType::N>(
67 nat_coords.data(), shape_matrix, GlobalDim,
68 false /* axisymmetric */);
69
70 matched = true;
71 }
72 });
73
74 if (!matched)
75 {
76 OGS_FATAL("Element type '{:s}' is not supported as anchor element.",
78 }
79
80 return shape_matrix.N;
81}
82
83template <int GlobalDim>
84std::tuple<Eigen::VectorXd, Eigen::MatrixXd> assembleLocalBJac(
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)
90{
91 // Signs for the two nodes alternate.
92 constexpr auto even_odd_sign = [](std::size_t const n)
93 { return (n % 2 == 0) ? 1.0 : -1.0; };
94
95 Eigen::VectorXd local_rhs(num_dof);
96 Eigen::MatrixXd local_Jac(num_dof, num_dof);
97 local_rhs.setZero();
98 local_Jac.setZero();
99
100 for (std::size_t node_idx = 0; node_idx < 2; ++node_idx)
101 {
102 // Select the correct block according to the node index.
103 // Replicate shape matrices of the corresponding element to GlobalDim
104 // dimensions. Multiply component-wise with the force vector f.
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]
109 .transpose()
110 .replicate<GlobalDim, 1>()
111 .cwiseProduct(f.transpose()
112 .replicate(nodes_per_element[node_idx], 1)
113 .reshaped());
114 for (std::size_t node_idx_inner = 0; node_idx_inner < 2;
115 ++node_idx_inner)
116 {
117 Eigen::MatrixXd const& ones = Eigen::MatrixXd::Ones(
118 nodes_per_element[node_idx], nodes_per_element[node_idx_inner]);
119
120 // Select the correct block according to the node index.
121 // Replicate shape matrices of the corresponding element to
122 // GlobalDim dimensions. Multiply component-wise with the derivative
123 // Df of force vector f.
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]
130 .transpose()
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));
135 }
136 }
137 return {local_rhs, local_Jac};
138}
139
140template <int GlobalDim>
142 MeshLib::Mesh const& bulk_mesh,
143 NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
144 std::size_t const source_term_mesh_id,
145 MeshLib::Mesh const& st_mesh,
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),
150 st_mesh_(st_mesh),
151 variable_id_(variable_id),
152 component_ids_(
153 []
154 {
155 std::array<int, GlobalDim> arr{};
156 std::iota(arr.begin(), arr.end(), 0);
157 return arr;
158 }())
159{
160 DBUG("Create EmbeddedAnchor.");
161 std::string_view const bulk_element_ids_string = "bulk_element_ids";
162 bulk_element_ids_ =
163 st_mesh_.getProperties().template getPropertyVector<std::size_t>(
164 bulk_element_ids_string);
165
166 std::string_view const natural_coordinates_string = "natural_coordinates";
167 natural_coordinates_ =
168 st_mesh_.getProperties().template getPropertyVector<double>(
169 natural_coordinates_string);
170
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);
176
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);
182
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);
188
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);
194
195 std::string_view const anchor_stiffness_string = "anchor_stiffness";
196 anchor_stiffness_ =
197 st_mesh_.getProperties().template getPropertyVector<double>(
198 anchor_stiffness_string);
199}
200
201template <int GlobalDim>
204 MeshLib::Element const* const anchor_element,
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,
209 GlobalVector const& x,
211{
212 for (auto&& [node_idx, anchor_node_id] : ranges::views::enumerate(
213 anchor_element->nodes() | MeshLib::views::ids))
214 {
215 // used pos in integrate function will be of the last node in anchor
216 // element, i,e, only one uniform stiffness is supported
217 auto const bulk_element_id = (*bulk_element_ids_)[anchor_node_id];
218 pos.setElementID(bulk_element_id);
219
220 auto const& bulk_element = *bulk_mesh_.getElement(bulk_element_id);
221 nodes_per_element[node_idx] = bulk_element.nodes().size();
222
223 std::span<double const, 3> const nat_coords(
224 &(*natural_coordinates_)[3 * anchor_node_id], 3);
225 auto N = computeShapeMatrix<GlobalDim>(bulk_element, nat_coords);
226 shape_matrices.push_back(std::move(N));
227
228 for (int component = 0; component < GlobalDim; ++component)
229 {
230 auto const& global_indices_percomponent = getIndices(
231 variable_id_, component, bulk_element_id, dof_table_bulk_);
232
233 Eigen::VectorXd const local_x_element_percomponent =
234 MathLib::toVector(x.get(global_indices_percomponent));
235
236 global_indices.insert(global_indices.cend(),
237 global_indices_percomponent.begin(),
238 global_indices_percomponent.end());
239
240 // this maps the displacement of each component of all nodes
241 // of the element to the desired anchor location in that element
242 // where all even local_x entries are for the first node and all odd
243 // local_x entries are for the second node
245 local_x_element_percomponent,
246 shape_matrices[node_idx],
247 local_x[nodeLocalIndices<GlobalDim>(node_idx)[component]]);
248 }
249 }
250}
251
252template <int GlobalDim>
254 GlobalVector const& x,
255 GlobalVector& b,
256 GlobalMatrix* jac) const
257{
258 DBUG("Assemble EmbeddedAnchor.");
259
260 using GlobalDimVector = Eigen::Vector<double, GlobalDim>;
261 using GlobalDimMatrix = Eigen::Matrix<double, GlobalDim, GlobalDim>;
262
263 for (MeshLib::Element const* const anchor_element : st_mesh_.getElements())
264 {
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,
273 local_x, x, pos);
274
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();
280
281 // Displacement in the two nodes.
282 auto u = [&local_x](int const i)
283 { return local_x(nodeLocalIndices<GlobalDim>(i)); };
284 GlobalDimVector const l = l_original + u(1) - u(0);
285
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];
296
297 double const strain = (l.norm() - l_original_norm) / l_original_norm;
298
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;
305
306 GlobalDimMatrix const Df_friction = GlobalDimMatrix::Zero();
307 GlobalDimMatrix const Df_elastic = l_original / l_original_norm * K *
308 l.transpose() / l.norm() /
309 l_original_norm;
310 GlobalDimMatrix const Df =
311 (f_elastic.norm() < max_force) ? Df_elastic : Df_friction;
312
313 auto const& [local_rhs, local_Jac] = assembleLocalBJac<GlobalDim>(
314 f, Df, shape_matrices, global_indices.size(), nodes_per_element);
315
316 b.add(global_indices, local_rhs);
317 if (jac)
318 {
319 jac->add({global_indices, global_indices}, local_Jac);
320 }
321 }
322}
323
324template class EmbeddedAnchor<2>;
325template class EmbeddedAnchor<3>;
326
327} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
int add(IndexType row, IndexType col, double val)
Definition EigenMatrix.h:87
Global vector based on Eigen vector.
Definition EigenVector.h:26
void add(IndexType rowId, double v)
add entry
Definition EigenVector.h:77
double get(IndexType rowId) const
get entry
Definition EigenVector.h:59
virtual CellType getCellType() const =0
virtual unsigned getNumberOfNodes() const =0
virtual constexpr unsigned getDimension() const =0
Get dimension of the mesh element.
constexpr std::span< Node *const > nodes() const
Span of element's nodes, their pointers actually.
Definition Element.h:74
void setElementID(std::size_t element_id)
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
void integrate(const double t, GlobalVector const &x, GlobalVector &b, GlobalMatrix *jac) const override
EmbeddedAnchor(MeshLib::Mesh const &bulk_mesh, NumLib::LocalToGlobalIndexMap const &dof_table_bulk, std::size_t const source_term_mesh_id, MeshLib::Mesh const &st_mesh, const int variable_id)
void foreach(Function &&f)
Definition TMP.h:157
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
Definition Mesh.h:227
std::string CellType2String(const CellType t)
Given a MeshElemType this returns the appropriate string.
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
NumLib::TemplateIsoparametric< ShapeFunction, ShapeMatricesType > createIsoparametricFiniteElement(MeshLib::Element const &e)
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)
Eigen::RowVectorXd computeShapeMatrix(MeshLib::Element const &bulk_element, std::span< double const, 3 > const nat_coords)
auto nodeLocalIndices(std::size_t const i)