OGS 6.1.0-1699-ge946d4c5f
SmallDeformationNonlocalProcess.cpp
Go to the documentation of this file.
1 
13 
14 #include <nlohmann/json.hpp>
15 #include <iostream>
16 
17 // Reusing local assembler creation code.
19 
20 namespace ProcessLib
21 {
22 namespace SmallDeformationNonlocal
23 {
24 template <int DisplacementDim>
27  MeshLib::Mesh& mesh,
28  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
29  jacobian_assembler,
30  std::vector<std::unique_ptr<ParameterBase>> const& parameters,
31  unsigned const integration_order,
32  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
33  process_variables,
35  SecondaryVariableCollection&& secondary_variables,
36  NumLib::NamedFunctionCaller&& named_function_caller)
37  : Process(mesh, std::move(jacobian_assembler), parameters,
38  integration_order, std::move(process_variables),
39  std::move(secondary_variables), std::move(named_function_caller)),
40  _process_data(std::move(process_data))
41 {
42  _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
43  mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
44 
45  _integration_point_writer.emplace_back(
46  std::make_unique<SigmaIntegrationPointWriter>(
47  static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
48  integration_order, [this]() {
49  // Result containing integration point data for each local
50  // assembler.
51  std::vector<std::vector<double>> result;
52  result.resize(_local_assemblers.size());
53 
54  for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
55  {
56  auto const& local_asm = *_local_assemblers[i];
57 
58  result[i] = local_asm.getSigma();
59  }
60 
61  return result;
62  }));
63 
64  _integration_point_writer.emplace_back(
65  std::make_unique<KappaDIntegrationPointWriter>(
66  integration_order, [this]() {
67  // Result containing integration point data for each local
68  // assembler.
69  std::vector<std::vector<double>> result;
70  result.resize(_local_assemblers.size());
71 
72  for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
73  {
74  auto const& local_asm = *_local_assemblers[i];
75 
76  result[i] = local_asm.getKappaD();
77  }
78 
79  return result;
80  }));
81 }
82 
83 template <int DisplacementDim>
85 {
86  return false;
87 }
88 
89 template <int DisplacementDim>
92  MeshLib::Mesh const& mesh,
93  unsigned const integration_order)
94 {
95  // Reusing local assembler creation code.
98  mesh.getElements(), dof_table, _local_assemblers,
99  mesh.isAxiallySymmetric(), integration_order, _process_data);
100 
101  // TODO move the two data members somewhere else.
102  // for extrapolation of secondary variables
103  std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
106  std::make_unique<NumLib::LocalToGlobalIndexMap>(
107  std::move(all_mesh_subsets_single_component),
108  // by location order is needed for output
110 
112  "sigma",
114  DisplacementDim>::RowsAtCompileTime,
115  getExtrapolator(), _local_assemblers,
117 
119  "epsilon",
121  DisplacementDim>::RowsAtCompileTime,
122  getExtrapolator(), _local_assemblers,
124 
126  "eps_p_V",
127  makeExtrapolator(1, getExtrapolator(), _local_assemblers,
130  "eps_p_D_xx",
131  makeExtrapolator(1, getExtrapolator(), _local_assemblers,
133 
135  "damage",
136  makeExtrapolator(1, getExtrapolator(), _local_assemblers,
138 
140  &LocalAssemblerInterface::nonlocal, _local_assemblers,
141  _local_assemblers);
142 
143  // Set initial conditions for integration point data.
144  for (auto const& ip_writer : _integration_point_writer)
145  {
146  auto const& name = ip_writer->name();
147  // First check the field data, which is used for restart.
148  if (mesh.getProperties().existsPropertyVector<double>(name))
149  {
150  auto const& mesh_property =
151  *mesh.getProperties().template getPropertyVector<double>(name);
152 
153  // The mesh property must be defined on integration points.
154  if (mesh_property.getMeshItemType() !=
156  {
157  continue;
158  }
159 
160  auto const ip_meta_data = getIntegrationPointMetaData(mesh, name);
161 
162  // Check the number of components.
163  if (ip_meta_data.n_components !=
164  mesh_property.getNumberOfComponents())
165  {
166  OGS_FATAL(
167  "Different number of components in meta data (%d) than in "
168  "the integration point field data for '%s': %d.",
169  ip_meta_data.n_components, name.c_str(),
170  mesh_property.getNumberOfComponents());
171  }
172 
173  // Now we have a properly named vtk's field data array and the
174  // corresponding meta data.
175  std::size_t position = 0;
176  for (auto& local_asm : _local_assemblers)
177  {
178  std::size_t const integration_points_read =
179  local_asm->setIPDataInitialConditions(
180  name, &mesh_property[position],
181  ip_meta_data.integration_order);
182  if (integration_points_read == 0)
183  {
184  OGS_FATAL(
185  "No integration points read in the integration point "
186  "initial conditions set function.");
187  }
188  position += integration_points_read * ip_meta_data.n_components;
189  }
190  }
191  else if (mesh.getProperties().existsPropertyVector<double>(name +
192  "_ic"))
193  { // Try to find cell data with '_ic' suffix
194  auto const& mesh_property =
195  *mesh.getProperties().template getPropertyVector<double>(name +
196  "_ic");
197  if (mesh_property.getMeshItemType() != MeshLib::MeshItemType::Cell)
198  {
199  continue;
200  }
201 
202  // Now we have a vtk's cell data array containing the initial
203  // conditions for the corresponding integration point writer.
204 
205  // For each assembler use the single cell value for all
206  // integration points.
207  for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
208  {
209  auto& local_asm = _local_assemblers[i];
210 
211  std::vector<double> value(
212  &mesh_property[i],
213  &mesh_property[i] + mesh_property.getNumberOfComponents());
214  // TODO (naumov) Check sizes / read size / etc.
215  // OR reconstruct dimensions from size / component =
216  // ip_points
217  local_asm->setIPDataInitialConditionsFromCellData(name, value);
218  }
219  }
220  }
221 }
222 
223 template <int DisplacementDim>
225  const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
226  GlobalVector& b)
227 {
228  DBUG("Assemble SmallDeformationNonlocalProcess.");
229 
230  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
231  dof_table = {std::ref(*_local_to_global_index_map)};
232 
233  const int process_id = 0;
234  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
235 
236  // Call global assembler for each local assembly item.
239  pv.getActiveElementIDs(), dof_table, t, x, M, K, b, _coupled_solutions);
240 }
241 
242 template <int DisplacementDim>
244  DisplacementDim>::preAssembleConcreteProcess(const double t,
245  GlobalVector const& x)
246 {
247  DBUG("preAssemble SmallDeformationNonlocalProcess.");
248 
249  const int process_id = 0;
250  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
251 
252  // Call global assembler for each local assembly item.
257 }
258 
259 template <int DisplacementDim>
261  assembleWithJacobianConcreteProcess(const double t, GlobalVector const& x,
262  GlobalVector const& xdot,
263  const double dxdot_dx,
264  const double dx_dx, GlobalMatrix& M,
265  GlobalMatrix& K, GlobalVector& b,
266  GlobalMatrix& Jac)
267 {
268  DBUG("AssembleWithJacobian SmallDeformationNonlocalProcess.");
269 
270  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
271  dof_table = {std::ref(*_local_to_global_index_map)};
272 
273  const int process_id = 0;
274  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
275 
276  // Call global assembler for each local assembly item.
279  _local_assemblers, pv.getActiveElementIDs(), dof_table, t, x, xdot,
280  dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
281 
282  b.copyValues(*_nodal_forces);
283  std::transform(_nodal_forces->begin(), _nodal_forces->end(),
284  _nodal_forces->begin(), [](double val) { return -val; });
285 }
286 
287 template <int DisplacementDim>
289  DisplacementDim>::preTimestepConcreteProcess(GlobalVector const& x,
290  double const t,
291  double const dt,
292  int const /*process_id*/)
293 {
294  DBUG("PreTimestep SmallDeformationNonlocalProcess.");
295 
296  _process_data.dt = dt;
297  _process_data.t = t;
298 
299  const int process_id = 0;
300  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
301 
305 }
306 
307 template <int DisplacementDim>
310  GlobalVector const& x)
311 {
312  _process_data.crack_volume_old = _process_data.crack_volume;
313  _process_data.crack_volume = 0.0;
314 
315  DBUG("PostNonLinearSolver crack volume computation.");
316 
317  const int process_id = 0;
318  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
319 
323  _process_data.crack_volume);
324 
325  INFO("Integral of crack: %g", _process_data.crack_volume);
326 
328 }
329 
332 
333 } // namespace SmallDeformationNonlocal
334 } // 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:262
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:123
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)
SmallDeformationNonlocalProcess(MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< 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)
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:131
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:150
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:134
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition: Process.h:264
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 preTimestep(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
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:266
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 std::vector< double > const & getIntPtEpsPV(const double, GlobalVector const &, NumLib::LocalToGlobalIndexMap const &, std::vector< double > &cache) const =0
CoupledSolutionsForStaggeredScheme * _coupled_solutions
Definition: Process.h:279
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:71
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:290
VectorMatrixAssembler _global_assembler
Definition: Process.h:273
Eigen::Matrix< double, KelvinVectorDimensions< DisplacementDim >::value, 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:49
std::vector< std::size_t > & getActiveElementIDs() const