OGS 6.2.1-97-g73d1aeda3
SmallDeformationProcess.cpp
Go to the documentation of this file.
1 
11 
12 #include <cassert>
13 #include <nlohmann/json.hpp>
14 
15 #include "BaseLib/Functional.h"
17 #include "ProcessLib/Process.h"
19 
20 #include "SmallDeformationFEM.h"
21 
22 namespace ProcessLib
23 {
24 namespace SmallDeformation
25 {
26 template <int DisplacementDim>
28  std::string name,
29  MeshLib::Mesh& mesh,
30  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
31  std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
32  unsigned const integration_order,
33  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
34  process_variables,
36  SecondaryVariableCollection&& secondary_variables,
37  NumLib::NamedFunctionCaller&& named_function_caller)
38  : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
39  integration_order, std::move(process_variables),
40  std::move(secondary_variables), std::move(named_function_caller)),
41  _process_data(std::move(process_data))
42 {
43  _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
44  mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
45 
46  _material_forces = MeshLib::getOrCreateMeshProperty<double>(
47  mesh, "MaterialForces", MeshLib::MeshItemType::Node, DisplacementDim);
48 
49  // TODO (naumov) remove ip suffix. Probably needs modification of the mesh
50  // properties, s.t. there is no "overlapping" with cell/point data.
51  // See getOrCreateMeshProperty.
52  _integration_point_writer.emplace_back(
53  std::make_unique<IntegrationPointWriter>(
54  "sigma_ip",
55  static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
56  integration_order, [this]() {
57  // Result containing integration point data for each local
58  // assembler.
59  std::vector<std::vector<double>> result;
60  result.resize(_local_assemblers.size());
61 
62  for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
63  {
64  auto const& local_asm = *_local_assemblers[i];
65 
66  result[i] = local_asm.getSigma();
67  }
68 
69  return result;
70  }));
71 }
72 
73 template <int DisplacementDim>
75 {
76  return false;
77 }
78 
79 template <int DisplacementDim>
81  NumLib::LocalToGlobalIndexMap const& dof_table,
82  MeshLib::Mesh const& mesh,
83  unsigned const integration_order)
84 {
85  using nlohmann::json;
86 
88  DisplacementDim, SmallDeformationLocalAssembler>(
89  mesh.getElements(), dof_table, _local_assemblers,
90  mesh.isAxiallySymmetric(), integration_order, _process_data);
91 
92  // TODO move the two data members somewhere else.
93  // for extrapolation of secondary variables
94  std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
97  std::make_unique<NumLib::LocalToGlobalIndexMap>(
98  std::move(all_mesh_subsets_single_component),
99  // by location order is needed for output
101 
103  "free_energy_density",
104  makeExtrapolator(1, getExtrapolator(), _local_assemblers,
106 
108  "sigma",
110  DisplacementDim>::RowsAtCompileTime,
111  getExtrapolator(), _local_assemblers,
113 
115  "epsilon",
117  DisplacementDim>::RowsAtCompileTime,
118  getExtrapolator(), _local_assemblers,
120 
121  //
122  // enable output of internal variables defined by material models
123  //
124 
125  // Collect the internal variables for all solid materials.
126  std::vector<typename MaterialLib::Solids::MechanicsBase<
127  DisplacementDim>::InternalVariable>
128  internal_variables;
129  for (auto const& material_id__solid_material :
130  _process_data.solid_materials)
131  {
132  auto const variables =
133  material_id__solid_material.second->getInternalVariables();
134  copy(begin(variables), end(variables),
135  back_inserter(internal_variables));
136  }
137 
138  // Register the internal variables.
139  for (auto const& internal_variable : internal_variables)
140  {
141  auto const& name = internal_variable.name;
142  auto const& fct = internal_variable.getter;
143  auto const num_components = internal_variable.num_components;
144  DBUG("Registering internal variable %s.", name.c_str());
145 
146  auto getIntPtValues = BaseLib::easyBind(
147  [fct, num_components](
148  LocalAssemblerInterface const& loc_asm,
149  const double /*t*/,
150  GlobalVector const& /*current_solution*/,
151  NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
152  std::vector<double>& cache) -> std::vector<double> const& {
153  const unsigned num_int_pts =
155 
156  cache.clear();
157  auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
158  double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(
159  cache, num_components, num_int_pts);
160 
161  // TODO avoid the heap allocation (one per finite element)
162  std::vector<double> cache_column(num_int_pts);
163 
164  for (unsigned i = 0; i < num_int_pts; ++i)
165  {
166  auto const& state = loc_asm.getMaterialStateVariablesAt(i);
167 
168  auto const& int_pt_values = fct(state, cache_column);
169  assert(int_pt_values.size() == num_components);
170  auto const int_pt_values_vec =
171  MathLib::toVector(int_pt_values);
172 
173  cache_mat.col(i).noalias() = int_pt_values_vec;
174  }
175 
176  return cache;
177  });
178 
180  name,
181  makeExtrapolator(num_components, getExtrapolator(),
182  _local_assemblers, std::move(getIntPtValues)));
183  }
184 
185  // Set initial conditions for integration point data.
186  for (auto const& ip_writer : _integration_point_writer)
187  {
188  // Find the mesh property with integration point writer's name.
189  auto const& name = ip_writer->name();
190  if (!mesh.getProperties().existsPropertyVector<double>(name))
191  {
192  continue;
193  }
194  auto const& mesh_property =
195  *mesh.getProperties().template getPropertyVector<double>(name);
196 
197  // The mesh property must be defined on integration points.
198  if (mesh_property.getMeshItemType() !=
200  {
201  continue;
202  }
203 
204  auto const ip_meta_data = getIntegrationPointMetaData(mesh, name);
205 
206  // Check the number of components.
207  if (ip_meta_data.n_components != mesh_property.getNumberOfComponents())
208  {
209  OGS_FATAL(
210  "Different number of components in meta data (%d) than in "
211  "the integration point field data for '%s': %d.",
212  ip_meta_data.n_components, name.c_str(),
213  mesh_property.getNumberOfComponents());
214  }
215 
216  // Now we have a properly named vtk's field data array and the
217  // corresponding meta data.
218  std::size_t position = 0;
219  for (auto& local_asm : _local_assemblers)
220  {
221  std::size_t const integration_points_read =
222  local_asm->setIPDataInitialConditions(
223  name, &mesh_property[position],
224  ip_meta_data.integration_order);
225  if (integration_points_read == 0)
226  {
227  OGS_FATAL(
228  "No integration points read in the integration point "
229  "initial conditions set function.");
230  }
231  position += integration_points_read * ip_meta_data.n_components;
232  }
233  }
234 
235  // Initialize local assemblers after all variables have been set.
237  &LocalAssemblerInterface::initialize, _local_assemblers,
239 }
240 
241 template <int DisplacementDim>
243  const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
244  GlobalVector& b)
245 {
246  DBUG("Assemble SmallDeformationProcess.");
247 
248  const int process_id = 0;
249  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
250 
251  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
252  dof_table = {std::ref(*_local_to_global_index_map)};
253  // Call global assembler for each local assembly item.
256  pv.getActiveElementIDs(), dof_table, t, x, M, K, b,
258 }
259 
260 template <int DisplacementDim>
262  assembleWithJacobianConcreteProcess(const double t, GlobalVector const& x,
263  GlobalVector const& xdot,
264  const double dxdot_dx,
265  const double dx_dx, GlobalMatrix& M,
266  GlobalMatrix& K, GlobalVector& b,
267  GlobalMatrix& Jac)
268 {
269  DBUG("AssembleWithJacobian SmallDeformationProcess.");
270 
271  const int process_id = 0;
272  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
273 
274  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
275  dof_table = {std::ref(*_local_to_global_index_map)};
276  // Call global assembler for each local assembly item.
279  _local_assemblers, pv.getActiveElementIDs(), dof_table, t, x,
280  xdot, dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
281 
283  *_nodal_forces, std::negate<double>());
284 }
285 
286 template <int DisplacementDim>
288  GlobalVector const& /*x*/, double const t, double const dt,
289  const int /*process_id*/)
290 {
291  DBUG("PreTimestep SmallDeformationProcess.");
292 
293  _process_data.dt = dt;
294  _process_data.t = t;
295 }
296 
297 template <int DisplacementDim>
299  GlobalVector const& x, const double /*t*/, const double /*delta_t*/,
300  int const process_id)
301 {
302  DBUG("PostTimestep SmallDeformationProcess.");
303 
304  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
305 
309 
310  std::unique_ptr<GlobalVector> material_forces;
312  material_forces, _local_assemblers, *_local_to_global_index_map, x);
313 
314  material_forces->copyValues(*_material_forces);
315 }
316 
317 template class SmallDeformationProcess<2>;
318 template class SmallDeformationProcess<3>;
319 
320 } // namespace SmallDeformation
321 } // namespace ProcessLib
IntegrationPointMetaData getIntegrationPointMetaData(MeshLib::Mesh const &mesh, std::string const &name)
virtual std::vector< InternalVariable > getInternalVariables() const
std::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_all_nodes
Definition: Process.h:288
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition: Process.h:125
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map_single_component
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
static void executeMemberOnDereferenced(Method method, Container const &container, Args &&... args)
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, double const t, GlobalVector const &x, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, CoupledSolutionsForStaggeredScheme const *const cpl_xs)
MeshLib::Properties & getProperties()
Definition: Mesh.h:134
void createLocalAssemblers(std::vector< MeshLib::Element *> const &mesh_elements, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< std::unique_ptr< LocalAssemblerInterface >> &local_assemblers, ExtraCtorArgs &&... extra_ctor_args)
NumLib::Extrapolator & getExtrapolator() const
Definition: Process.h:152
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, const double t, GlobalVector const &x, GlobalVector const &xdot, const double dxdot_dx, const double dx_dx, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac, CoupledSolutionsForStaggeredScheme const *const cpl_xs)
bool isAxiallySymmetric() const
Definition: Mesh.h:137
void preTimestepConcreteProcess(GlobalVector const &x, double const t, double const dt, const int process_id) override
virtual std::vector< double > const & getIntPtFreeEnergyDensity(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const =0
virtual void initialize(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
SmallDeformationProcess(std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase >> const &parameters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable >>> &&process_variables, SmallDeformationProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, NumLib::NamedFunctionCaller &&named_function_caller)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:31
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition: Process.h:290
Builds expression trees of named functions dynamically at runtime.
virtual std::vector< double > const & getIntPtEpsilon(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const =0
void assembleConcreteProcess(const double t, GlobalVector const &x, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
SecondaryVariableCollection _secondary_variables
Definition: Process.h:292
std::vector< std::unique_ptr< LocalAssemblerInterface > > _local_assemblers
void writeMaterialForces(std::unique_ptr< GlobalVector > &material_forces, std::vector< std::unique_ptr< LocalAssemblerInterface >> const &local_assemblers, NumLib::LocalToGlobalIndexMap const &local_to_global_index_map, GlobalVector const &x)
std::enable_if< std::is_same< MethodClass, typename std::remove_cv< typename std::remove_pointer< typename std::decay< Object >::type >::type >::type >::value, std::function< ReturnType(Args...)> >::type easyBind(ReturnType(MethodClass::*method)(Args...), Object &&obj)
Definition: Functional.h:180
virtual MaterialLib::Solids::MechanicsBase< DisplacementDim >::MaterialStateVariables const & getMaterialStateVariablesAt(unsigned) const =0
void postTimestepConcreteProcess(GlobalVector const &x, const double t, const double delta_t, int const process_id) override
SmallDeformationProcessData< DisplacementDim > _process_data
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:108
bool existsPropertyVector(std::string const &name) const
Definition: Properties.h:79
virtual void postTimestep(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x)
virtual std::vector< double > const & getIntPtSigma(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const =0
void transformVariableFromGlobalVector(GlobalVector const &input_vector, int const variable_id, NumLib::LocalToGlobalIndexMap const &local_to_global_index_map, MeshLib::PropertyVector< double > &output_vector, Functor mapFunction)
Definition: DOFTableUtil.h:59
CoupledSolutionsForStaggeredScheme * _coupled_solutions
Definition: Process.h:305
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
Handles configuration of several secondary variables from the project file.
Ordering data by spatial location.
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
void addSecondaryVariable(std::string const &internal_name, SecondaryVariableFunctions &&fcts)
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:36
void assembleWithJacobianConcreteProcess(const double t, GlobalVector const &x, GlobalVector const &xdot, const double dxdot_dx, const double dx_dx, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
std::vector< std::unique_ptr< IntegrationPointWriter > > _integration_point_writer
Definition: Process.h:316
VectorMatrixAssembler _global_assembler
Definition: Process.h:299
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
std::string const name
Definition: Process.h:284
std::vector< std::size_t > & getActiveElementIDs() const