OGS 6.2.1-97-g73d1aeda3
SmallDeformationNonlocalProcess.cpp
Go to the documentation of this file.
1 
13 
14 #include <nlohmann/json.hpp>
15 #include <iostream>
16 
18 
19 // Reusing local assembler creation code.
21 
22 namespace ProcessLib
23 {
24 namespace SmallDeformationNonlocal
25 {
26 template <int DisplacementDim>
29  std::string name,
30  MeshLib::Mesh& mesh,
31  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
32  jacobian_assembler,
33  std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
34  parameters,
35  unsigned const integration_order,
36  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
37  process_variables,
39  SecondaryVariableCollection&& secondary_variables,
40  NumLib::NamedFunctionCaller&& named_function_caller)
41  : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
42  integration_order, std::move(process_variables),
43  std::move(secondary_variables), std::move(named_function_caller)),
44  _process_data(std::move(process_data))
45 {
46  _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
47  mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
48 
49  _integration_point_writer.emplace_back(
50  std::make_unique<IntegrationPointWriter>(
51  "sigma_ip",
52  static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
53  integration_order, [this]() {
54  // Result containing integration point data for each local
55  // assembler.
56  std::vector<std::vector<double>> result;
57  result.resize(_local_assemblers.size());
58 
59  for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
60  {
61  auto const& local_asm = *_local_assemblers[i];
62 
63  result[i] = local_asm.getSigma();
64  }
65 
66  return result;
67  }));
68 
69  _integration_point_writer.emplace_back(
70  std::make_unique<IntegrationPointWriter>(
71  "kappa_d_ip", 1 /*n_components*/, integration_order, [this]() {
72  // Result containing integration point data for each local
73  // assembler.
74  std::vector<std::vector<double>> result;
75  result.resize(_local_assemblers.size());
76 
77  for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
78  {
79  auto const& local_asm = *_local_assemblers[i];
80 
81  result[i] = local_asm.getKappaD();
82  }
83 
84  return result;
85  }));
86 }
87 
88 template <int DisplacementDim>
90 {
91  return false;
92 }
93 
94 template <int DisplacementDim>
97  MeshLib::Mesh const& mesh,
98  unsigned const integration_order)
99 {
100  // Reusing local assembler creation code.
102  DisplacementDim, SmallDeformationNonlocalLocalAssembler>(
103  mesh.getElements(), dof_table, _local_assemblers,
104  mesh.isAxiallySymmetric(), integration_order, _process_data);
105 
106  // TODO move the two data members somewhere else.
107  // for extrapolation of secondary variables
108  std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
111  std::make_unique<NumLib::LocalToGlobalIndexMap>(
112  std::move(all_mesh_subsets_single_component),
113  // by location order is needed for output
115 
117  "sigma",
119  DisplacementDim>::RowsAtCompileTime,
120  getExtrapolator(), _local_assemblers,
122 
124  "epsilon",
126  DisplacementDim>::RowsAtCompileTime,
127  getExtrapolator(), _local_assemblers,
129 
131  "eps_p_V",
132  makeExtrapolator(1, getExtrapolator(), _local_assemblers,
135  "eps_p_D_xx",
136  makeExtrapolator(1, getExtrapolator(), _local_assemblers,
138 
140  "damage",
141  makeExtrapolator(1, getExtrapolator(), _local_assemblers,
143 
145  &LocalAssemblerInterface::nonlocal, _local_assemblers,
146  _local_assemblers);
147 
148  // Set initial conditions for integration point data.
149  for (auto const& ip_writer : _integration_point_writer)
150  {
151  auto const& name = ip_writer->name();
152  // First check the field data, which is used for restart.
153  if (mesh.getProperties().existsPropertyVector<double>(name))
154  {
155  auto const& mesh_property =
156  *mesh.getProperties().template getPropertyVector<double>(name);
157 
158  // The mesh property must be defined on integration points.
159  if (mesh_property.getMeshItemType() !=
161  {
162  continue;
163  }
164 
165  auto const ip_meta_data = getIntegrationPointMetaData(mesh, name);
166 
167  // Check the number of components.
168  if (ip_meta_data.n_components !=
169  mesh_property.getNumberOfComponents())
170  {
171  OGS_FATAL(
172  "Different number of components in meta data (%d) than in "
173  "the integration point field data for '%s': %d.",
174  ip_meta_data.n_components, name.c_str(),
175  mesh_property.getNumberOfComponents());
176  }
177 
178  // Now we have a properly named vtk's field data array and the
179  // corresponding meta data.
180  std::size_t position = 0;
181  for (auto& local_asm : _local_assemblers)
182  {
183  std::size_t const integration_points_read =
184  local_asm->setIPDataInitialConditions(
185  name, &mesh_property[position],
186  ip_meta_data.integration_order);
187  if (integration_points_read == 0)
188  {
189  OGS_FATAL(
190  "No integration points read in the integration point "
191  "initial conditions set function.");
192  }
193  position += integration_points_read * ip_meta_data.n_components;
194  }
195  }
196  else if (mesh.getProperties().existsPropertyVector<double>(name +
197  "_ic"))
198  { // Try to find cell data with '_ic' suffix
199  auto const& mesh_property =
200  *mesh.getProperties().template getPropertyVector<double>(name +
201  "_ic");
202  if (mesh_property.getMeshItemType() != MeshLib::MeshItemType::Cell)
203  {
204  continue;
205  }
206 
207  // Now we have a vtk's cell data array containing the initial
208  // conditions for the corresponding integration point writer.
209 
210  // For each assembler use the single cell value for all
211  // integration points.
212  for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
213  {
214  auto& local_asm = _local_assemblers[i];
215 
216  std::vector<double> value(
217  &mesh_property[i],
218  &mesh_property[i] + mesh_property.getNumberOfComponents());
219  // TODO (naumov) Check sizes / read size / etc.
220  // OR reconstruct dimensions from size / component =
221  // ip_points
222  local_asm->setIPDataInitialConditionsFromCellData(name, value);
223  }
224  }
225  }
226 
227  // Initialize local assemblers after all variables have been set.
229  &LocalAssemblerInterface::initialize, _local_assemblers,
231 }
232 
233 template <int DisplacementDim>
235  const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
236  GlobalVector& b)
237 {
238  DBUG("Assemble SmallDeformationNonlocalProcess.");
239 
240  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
241  dof_table = {std::ref(*_local_to_global_index_map)};
242 
243  const int process_id = 0;
244  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
245 
246  // Call global assembler for each local assembly item.
249  pv.getActiveElementIDs(), dof_table, t, x, M, K, b, _coupled_solutions);
250 }
251 
252 template <int DisplacementDim>
254  DisplacementDim>::preAssembleConcreteProcess(const double t,
255  GlobalVector const& x)
256 {
257  DBUG("preAssemble SmallDeformationNonlocalProcess.");
258 
259  const int process_id = 0;
260  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
261 
262  // Call global assembler for each local assembly item.
267 }
268 
269 template <int DisplacementDim>
271  assembleWithJacobianConcreteProcess(const double t, GlobalVector const& x,
272  GlobalVector const& xdot,
273  const double dxdot_dx,
274  const double dx_dx, GlobalMatrix& M,
275  GlobalMatrix& K, GlobalVector& b,
276  GlobalMatrix& Jac)
277 {
278  DBUG("AssembleWithJacobian SmallDeformationNonlocalProcess.");
279 
280  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
281  dof_table = {std::ref(*_local_to_global_index_map)};
282 
283  const int process_id = 0;
284  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
285 
286  // Call global assembler for each local assembly item.
289  _local_assemblers, pv.getActiveElementIDs(), dof_table, t, x, xdot,
290  dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
291 
292  b.copyValues(*_nodal_forces);
293  std::transform(_nodal_forces->begin(), _nodal_forces->end(),
294  _nodal_forces->begin(), [](double val) { return -val; });
295 }
296 
297 template <int DisplacementDim>
299  DisplacementDim>::preTimestepConcreteProcess(GlobalVector const& /*x*/,
300  double const t,
301  double const dt,
302  int const /*process_id*/)
303 {
304  DBUG("PreTimestep SmallDeformationNonlocalProcess.");
305 
306  _process_data.dt = dt;
307  _process_data.t = t;
308 }
309 
310 template <int DisplacementDim>
312  DisplacementDim>::postTimestepConcreteProcess(GlobalVector const& x,
313  double const /*t*/,
314  double const /*dt*/,
315  int const process_id)
316 {
317  DBUG("PostTimestep SmallDeformationNonlocalProcess.");
318 
319  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
320 
324 }
325 
326 template <int DisplacementDim>
329  GlobalVector const& x)
330 {
331  _process_data.crack_volume_old = _process_data.crack_volume;
332  _process_data.crack_volume = 0.0;
333 
334  DBUG("PostNonLinearSolver crack volume computation.");
335 
336  const int process_id = 0;
337  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
338 
342  _process_data.crack_volume);
343 
344  INFO("Integral of crack: %g", _process_data.crack_volume);
345 
347 }
348 
351 
352 } // namespace SmallDeformationNonlocal
353 } // namespace ProcessLib
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
IntegrationPointMetaData getIntegrationPointMetaData(MeshLib::Mesh const &mesh, std::string const &name)
SmallDeformationNonlocalProcessData< DisplacementDim > _process_data
virtual std::vector< double > const & getIntPtEpsPDXX(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const =0
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)
void assembleConcreteProcess(const double t, GlobalVector const &x, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
void preAssemble(const std::size_t mesh_item_id, LocalAssemblerInterface &local_assembler, const NumLib::LocalToGlobalIndexMap &dof_table, const double t, const GlobalVector &x)
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map_single_component
NumLib::IterationResult postIterationConcreteProcess(GlobalVector const &x) override
virtual std::vector< double > const & getIntPtEpsilon(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const =0
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)
void preTimestepConcreteProcess(GlobalVector const &x, double const t, double const dt, int const) override
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
virtual void initialize(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
SmallDeformationNonlocalProcess(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, SmallDeformationNonlocalProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, NumLib::NamedFunctionCaller &&named_function_caller)
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition: Process.h:290
IterationResult
Status flags telling the NonlinearSolver if an iteration succeeded.
virtual std::vector< double > const & getIntPtDamage(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const =0
Builds expression trees of named functions dynamically at runtime.
virtual void computeCrackIntegral(std::size_t mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double &crack_volume)=0
void preAssembleConcreteProcess(const double t, GlobalVector const &x) 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
void postTimestepConcreteProcess(GlobalVector const &x, double const t, double const dt, int const process_id) override
std::vector< std::unique_ptr< LocalAssemblerInterface > > _local_assemblers
virtual std::vector< double > const & getIntPtSigma(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const =0
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 & getIntPtEpsPV(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const =0
CoupledSolutionsForStaggeredScheme * _coupled_solutions
Definition: Process.h:305
virtual void nonlocal(std::size_t const mesh_item_id, std::vector< std::unique_ptr< SmallDeformationNonlocalLocalAssemblerInterface >> const &local_assemblers)=0
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 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