Merge branch 'anchor_elements' into 'master'

[PL/ST] implement anchor elements

See merge request ogs/ogs!5094

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Anchor/c_Anchor.md

@@ -0,0 +1,3 @@
+This source term models anchors that begin and end in a node of the bulk mesh.
+
+Each one-dimensional line element in the `Anchor` source term mesh defines an anchor.

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Anchor/t_anchor_force_constant.md

@@ -0,0 +1 @@
+Defines the parameter name for providing the anchor force constant (SI unit: `N`) between two nodes of an anchor element.

ProcessLib/BoundaryConditionAndSourceTerm/AnchorTerm.cpp

@@ -0,0 +1,118 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "AnchorTerm.h"
+
+#include <Eigen/Core>
+#include <cassert>
+
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+
+namespace ProcessLib
+{
+template <int GlobalDim>
+AnchorTerm<GlobalDim>::AnchorTerm(
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap> source_term_dof_table,
+ std::size_t const source_term_mesh_id,
+ MeshLib::Mesh const& st_mesh,
+ const int variable_id,
+ ParameterLib::Parameter<double> const& parameter)
+ : SourceTerm(std::move(source_term_dof_table)),
+ source_term_mesh_id_(source_term_mesh_id),
+ st_mesh_(st_mesh),
+ variable_id_(variable_id),
+ parameter_(parameter)
+{
+ DBUG("Create AnchorTerm.");
+}
+
+template <int GlobalDim>
+void AnchorTerm<GlobalDim>::integrate(const double t, GlobalVector const& x,
+ GlobalVector& b, GlobalMatrix* jac) const
+{
+ DBUG("Assemble AnchorTerm.");
+
+ using GlobalDimVector = Eigen::Vector<double, GlobalDim>;
+ using GlobalDimMatrix = Eigen::Matrix<double, GlobalDim, GlobalDim>;
+
+ for (MeshLib::Element const* const element : st_mesh_.getElements())
+ {
+ auto const element_id = element->getID();
+
+ ParameterLib::SpatialPosition pos;
+ pos.setElementID(element_id);
+
+ std::vector<GlobalIndexType> const global_indices =
+ NumLib::getIndices(element_id, *_source_term_dof_table);
+ assert(global_indices.size() == 2 * GlobalDim);
+
+ Eigen::Vector<double, 2 * GlobalDim> const local_x =
+ MathLib::toVector(x.get(global_indices));
+
+ Eigen::Vector<double, 2 * GlobalDim> local_rhs =
+ Eigen::Vector<double, 2 * GlobalDim>::Zero();
+ Eigen::Matrix<double, 2 * GlobalDim, 2 * GlobalDim> local_Jac =
+ Eigen::Matrix<double, 2 * GlobalDim, 2 * GlobalDim>::Zero();
+
+ // The local indices are, due to the nature of the DOF table, all even
+ // for the first node and odd for the second node.
+ auto node_local_indices = [](int const i)
+ { return Eigen::seqN(i, Eigen::fix<GlobalDim>, Eigen::fix<2>); };
+
+ auto node_coords = [element](int const i)
+ { return element->getNode(i)->asEigenVector3d(); };
+ GlobalDimVector const l_original =
+ (node_coords(1) - node_coords(0)).template head<GlobalDim>();
+ double const l_original_norm = l_original.norm();
+
+ // Displacement in the two nodes.
+ auto u = [&local_x, &node_local_indices](int const i)
+ { return local_x(node_local_indices(i)); };
+ GlobalDimVector const l = l_original + u(1) - u(0);
+
+ double const K = parameter_(t, pos)[0];
+ GlobalDimVector const f = l_original / l_original_norm * K *
+ (l.norm() - l_original_norm) /
+ l_original_norm;
+
+ GlobalDimMatrix const Df = l_original / l_original_norm * K *
+ l.transpose() / l.norm() / l_original_norm;
+
+ // Signs for the two nodes alternate.
+ constexpr auto even_odd_sign = [](int const n)
+ { return (n % 2 == 0) ? 1.0 : -1.0; };
+
+ // There are two blocks in the rhs vector for the two nodes of the
+ // element and four blocks in the local Jacobian formed by the four
+ // combinations of the two nodes.
+ for (int i = 0; i < 2; ++i)
+ {
+ local_rhs(node_local_indices(i)).noalias() += even_odd_sign(i) * f;
+
+ for (int j = 0; j < 2; ++j)
+ {
+ local_Jac(node_local_indices(i), node_local_indices(j))
+ .noalias() += even_odd_sign(i) * even_odd_sign(j) * Df;
+ }
+ }
+
+ b.add(global_indices, local_rhs);
+ if (jac)
+ {
+ jac->add({global_indices, global_indices}, local_Jac);
+ }
+ }
+}
+
+template class AnchorTerm<2>;
+template class AnchorTerm<3>;
+
+} // namespace ProcessLib

ProcessLib/BoundaryConditionAndSourceTerm/AnchorTerm.h

@@ -0,0 +1,39 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "SourceTerm.h"
+
+namespace ProcessLib
+{
+template <int GlobalDim>
+class AnchorTerm final : public SourceTerm
+{
+public:
+ explicit AnchorTerm(
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap> source_term_dof_table,
+ std::size_t const source_term_mesh_id, MeshLib::Mesh const& st_mesh,
+ const int variable_id,
+ ParameterLib::Parameter<double> const& parameter);
+
+ void integrate(const double t, GlobalVector const& x, GlobalVector& b,
+ GlobalMatrix* jac) const override;
+
+private:
+ std::size_t const source_term_mesh_id_;
+ MeshLib::Mesh const& st_mesh_;
+ int const variable_id_;
+ ParameterLib::Parameter<double> const& parameter_;
+};
+
+extern template class AnchorTerm<2>;
+extern template class AnchorTerm<3>;
+} // namespace ProcessLib

ProcessLib/BoundaryConditionAndSourceTerm/CreateAnchorTerm.cpp

@@ -0,0 +1,65 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreateAnchorTerm.h"
+
+#include "AnchorTerm.h"
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/Logging.h"
+#include "ParameterLib/Utils.h"
+
+namespace ProcessLib
+{
+template <int GlobalDim>
+std::unique_ptr<SourceTerm> createAnchorTerm(
+ BaseLib::ConfigTree const& config, MeshLib::Mesh const& st_mesh,
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap> dof_table,
+ std::size_t const source_term_mesh_id, const int variable_id,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters)
+{
+ DBUG("Constructing AnchorTerm from config.");
+ //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__type}
+ config.checkConfigParameter("type", "Anchor");
+
+ auto const param_name =
+ //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__Anchor__anchor_force_constant}
+ config.getConfigParameter<std::string>("anchor_force_constant");
+ DBUG("Using parameter {:s} as anchor stiffness constant.", param_name);
+
+ auto& param = ParameterLib::findParameter<double>(param_name, parameters, 1,
+ &st_mesh);
+
+ for (MeshLib::Element const* const element : st_mesh.getElements())
+ {
+ if (element->getNumberOfNodes() != 2)
+ {
+ OGS_FATAL(
+ "Every anchor element needs to have precisely two nodes.");
+ }
+ }
+
+ return std::make_unique<AnchorTerm<GlobalDim>>(
+ std::move(dof_table), source_term_mesh_id, st_mesh, variable_id, param);
+}
+
+template std::unique_ptr<SourceTerm> createAnchorTerm<2>(
+ BaseLib::ConfigTree const& config, MeshLib::Mesh const& st_mesh,
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap> dof_table,
+ std::size_t const source_term_mesh_id, const int variable_id,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
+ parameters);
+
+template std::unique_ptr<SourceTerm> createAnchorTerm<3>(
+ BaseLib::ConfigTree const& config, MeshLib::Mesh const& st_mesh,
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap> dof_table,
+ std::size_t const source_term_mesh_id, const int variable_id,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
+ parameters);
+} // namespace ProcessLib

ProcessLib/BoundaryConditionAndSourceTerm/CreateAnchorTerm.h

@@ -0,0 +1,47 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+namespace MeshLib
+{
+class Mesh;
+}
+namespace NumLib
+{
+class LocalToGlobalIndexMap;
+}
+namespace ParameterLib
+{
+struct ParameterBase;
+}
+namespace ProcessLib
+{
+class SourceTerm;
+}
+
+namespace ProcessLib
+{
+template <int GlobalDim>
+std::unique_ptr<SourceTerm> createAnchorTerm(
+ BaseLib::ConfigTree const& config, MeshLib::Mesh const& st_mesh,
+ std::unique_ptr<NumLib::LocalToGlobalIndexMap> dof_table,
+ std::size_t source_term_mesh_id, const int variable_id,
+ std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
+ parameters);
+
+} // namespace ProcessLib

ProcessLib/BoundaryConditionAndSourceTerm/CreateSourceTerm.cpp

@@ -10,6 +10,10 @@
 
 #include "CreateSourceTerm.h"
 
+#include <cassert>
+#include <numeric>
+
+#include "CreateAnchorTerm.h"
 #include "CreateNodalSourceTerm.h"
 #include "CreateVolumetricSourceTerm.h"
 #include "Python/CreatePythonSourceTerm.h"
@@ -75,7 +79,46 @@ std::unique_ptr<SourceTerm> createSourceTerm(
  source_term_mesh.getID(), variable_id, config.component_id,
  parameters);
  }
-
+ if (type == "Anchor")
+ {
+ const int number_of_components =
+ dof_table_bulk.getNumberOfVariableComponents(variable_id);
+ std::vector<int> component_ids(number_of_components);
+ std::iota(std::begin(component_ids), std::end(component_ids), 0);
+ auto dof_table_source_term =
+ dof_table_bulk.deriveBoundaryConstrainedMap(
+ variable_id, component_ids, std::move(source_term_mesh_subset));
+ const int bulk_mesh_dimension =
+ dof_table_bulk.getMeshSubset(variable_id, 0)
+ .getMesh()
+ .getDimension();
+ if (bulk_mesh_dimension != number_of_components)
+ {
+ OGS_FATAL(
+ "For the Anchor source term type,"
+ "the bulk mesh dimension needs to be the same "
+ "as the number of process variable components.");
+ }
+ switch (bulk_mesh_dimension)
+ {
+ case 2:
+ return ProcessLib::createAnchorTerm<2>(
+ config.config, config.mesh,
+ std::move(dof_table_source_term), source_term_mesh.getID(),
+ variable_id, parameters);
+ case 3:
+ return ProcessLib::createAnchorTerm<3>(
+ config.config, config.mesh,
+ std::move(dof_table_source_term), source_term_mesh.getID(),
+ variable_id, parameters);
+ default:
+ OGS_FATAL(
+ "Anchor can not be instantiated "
+ "for mesh dimensions other than two or three. "
+ "{}-dimensional mesh was given.",
+ bulk_mesh_dimension);
+ }
+ }
  if (type == "Line" || type == "Volumetric")
  {
  auto dof_table_source_term =
@@ -91,7 +134,6 @@ std::unique_ptr<SourceTerm> createSourceTerm(
  std::move(dof_table_source_term), parameters, integration_order,
  shapefunction_order);
  }
-
  if (type == "Python")
  {
  auto dof_table_source_term =