OGS
SmallDeformationNonlocalProcess.cpp
Go to the documentation of this file.
1 
12 
13 #include <iostream>
14 #include <nlohmann/json.hpp>
15 
17 
18 // Reusing local assembler creation code.
20 
21 namespace ProcessLib
22 {
23 namespace SmallDeformationNonlocal
24 {
25 template <int DisplacementDim>
28  std::string name,
29  MeshLib::Mesh& mesh,
30  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&&
31  jacobian_assembler,
32  std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
33  parameters,
34  unsigned const integration_order,
35  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
36  process_variables,
38  SecondaryVariableCollection&& secondary_variables)
39  : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
40  integration_order, std::move(process_variables),
41  std::move(secondary_variables)),
42  _process_data(std::move(process_data))
43 {
44  _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
45  mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
46 
47  _integration_point_writer.emplace_back(
48  std::make_unique<IntegrationPointWriter>(
49  "sigma_ip",
50  static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
51  integration_order, _local_assemblers,
53 
54  _integration_point_writer.emplace_back(
55  std::make_unique<IntegrationPointWriter>(
56  "kappa_d_ip", 1 /*n_components*/, integration_order,
58 }
59 
60 template <int DisplacementDim>
62 {
63  return false;
64 }
65 
66 template <int DisplacementDim>
69  MeshLib::Mesh const& mesh,
70  unsigned const integration_order)
71 {
72  // Reusing local assembler creation code.
75  mesh.getElements(), dof_table, _local_assemblers,
76  NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
77  _process_data);
78 
79  // TODO move the two data members somewhere else.
80  // for extrapolation of secondary variables
81  std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
82  *_mesh_subset_all_nodes};
83  _local_to_global_index_map_single_component =
84  std::make_unique<NumLib::LocalToGlobalIndexMap>(
85  std::move(all_mesh_subsets_single_component),
86  // by location order is needed for output
88 
90  "sigma",
92  DisplacementDim>::RowsAtCompileTime,
93  getExtrapolator(), _local_assemblers,
94  &LocalAssemblerInterface::getIntPtSigma));
95 
97  "epsilon",
99  DisplacementDim>::RowsAtCompileTime,
100  getExtrapolator(), _local_assemblers,
101  &LocalAssemblerInterface::getIntPtEpsilon));
102 
104  "eps_p_V",
105  makeExtrapolator(1, getExtrapolator(), _local_assemblers,
106  &LocalAssemblerInterface::getIntPtEpsPV));
108  "eps_p_D_xx",
109  makeExtrapolator(1, getExtrapolator(), _local_assemblers,
110  &LocalAssemblerInterface::getIntPtEpsPDXX));
111 
113  "damage",
114  makeExtrapolator(1, getExtrapolator(), _local_assemblers,
115  &LocalAssemblerInterface::getIntPtDamage));
116 
118  &LocalAssemblerInterface::nonlocal, _local_assemblers,
119  _local_assemblers);
120 
121  // Set initial conditions for integration point data.
122  for (auto const& ip_writer : _integration_point_writer)
123  {
124  auto const& name = ip_writer->name();
125  // First check the field data, which is used for restart.
126  if (mesh.getProperties().existsPropertyVector<double>(name))
127  {
128  auto const& mesh_property =
129  *mesh.getProperties().template getPropertyVector<double>(name);
130 
131  // The mesh property must be defined on integration points.
132  if (mesh_property.getMeshItemType() !=
134  {
135  continue;
136  }
137 
138  auto const ip_meta_data =
140 
141  // Check the number of components.
142  if (ip_meta_data.n_components !=
143  mesh_property.getNumberOfGlobalComponents())
144  {
145  OGS_FATAL(
146  "Different number of components in meta data ({:d}) than "
147  "in the integration point field data for '{:s}': {:d}.",
148  ip_meta_data.n_components, name,
149  mesh_property.getNumberOfGlobalComponents());
150  }
151 
152  // Now we have a properly named vtk's field data array and the
153  // corresponding meta data.
154  std::size_t position = 0;
155  for (auto& local_asm : _local_assemblers)
156  {
157  std::size_t const integration_points_read =
158  local_asm->setIPDataInitialConditions(
159  name, &mesh_property[position],
160  ip_meta_data.integration_order);
161  if (integration_points_read == 0)
162  {
163  OGS_FATAL(
164  "No integration points read in the integration point "
165  "initial conditions set function.");
166  }
167  position += integration_points_read * ip_meta_data.n_components;
168  }
169  }
170  else if (mesh.getProperties().existsPropertyVector<double>(name +
171  "_ic"))
172  { // Try to find cell data with '_ic' suffix
173  auto const& mesh_property =
174  *mesh.getProperties().template getPropertyVector<double>(name +
175  "_ic");
176  if (mesh_property.getMeshItemType() != MeshLib::MeshItemType::Cell)
177  {
178  continue;
179  }
180 
181  // Now we have a vtk's cell data array containing the initial
182  // conditions for the corresponding integration point writer.
183 
184  // For each assembler use the single cell value for all
185  // integration points.
186  for (std::size_t i = 0; i < _local_assemblers.size(); ++i)
187  {
188  auto& local_asm = _local_assemblers[i];
189 
190  std::vector<double> value(
191  &mesh_property[i],
192  &mesh_property[i] +
193  mesh_property.getNumberOfGlobalComponents());
194  // TODO (naumov) Check sizes / read size / etc.
195  // OR reconstruct dimensions from size / component =
196  // ip_points
197  local_asm->setIPDataInitialConditionsFromCellData(name, value);
198  }
199  }
200  }
201 
202  // Initialize local assemblers after all variables have been set.
204  &LocalAssemblerInterface::initialize, _local_assemblers,
205  *_local_to_global_index_map);
206 }
207 
208 template <int DisplacementDim>
210  const double t, double const dt, std::vector<GlobalVector*> const& x,
211  std::vector<GlobalVector*> const& xdot, int const process_id,
213 {
214  DBUG("Assemble SmallDeformationNonlocalProcess.");
215 
216  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
217  dof_table = {std::ref(*_local_to_global_index_map)};
218 
219  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
220 
221  // Call global assembler for each local assembly item.
223  _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
224  pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
225  b);
226 }
227 
228 template <int DisplacementDim>
230  DisplacementDim>::preAssembleConcreteProcess(const double t,
231  double const dt,
232  GlobalVector const& x)
233 {
234  DBUG("preAssemble SmallDeformationNonlocalProcess.");
235 
236  const int process_id = 0;
237  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
238 
239  // Call global assembler for each local assembly item.
241  _global_assembler, &VectorMatrixAssembler::preAssemble,
242  _local_assemblers, pv.getActiveElementIDs(),
243  *_local_to_global_index_map, t, dt, x);
244 }
245 
246 template <int DisplacementDim>
248  assembleWithJacobianConcreteProcess(const double t, double const dt,
249  std::vector<GlobalVector*> const& x,
250  std::vector<GlobalVector*> const& xdot,
251  int const process_id, GlobalMatrix& M,
252  GlobalMatrix& K, GlobalVector& b,
253  GlobalMatrix& Jac)
254 {
255  DBUG("AssembleWithJacobian SmallDeformationNonlocalProcess.");
256 
257  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
258  dof_table = {std::ref(*_local_to_global_index_map)};
259 
260  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
261 
262  // Call global assembler for each local assembly item.
265  _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
266  process_id, M, K, b, Jac);
267 
268  b.copyValues(*_nodal_forces);
269  std::transform(_nodal_forces->begin(), _nodal_forces->end(),
270  _nodal_forces->begin(), [](double val) { return -val; });
271 }
272 
273 template <int DisplacementDim>
275  postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
276  double const t,
277  double const dt,
278  int const process_id)
279 {
280  DBUG("PostTimestep SmallDeformationNonlocalProcess.");
281  std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
282  dof_tables.reserve(x.size());
283  std::generate_n(std::back_inserter(dof_tables), x.size(),
284  [&]() { return _local_to_global_index_map.get(); });
285 
286  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
288  &LocalAssemblerInterface::postTimestep, _local_assemblers,
289  pv.getActiveElementIDs(), dof_tables, x, t, dt);
290 }
291 
292 template <int DisplacementDim>
295  GlobalVector const& x)
296 {
297  _process_data.crack_volume_old = _process_data.crack_volume;
298  _process_data.crack_volume = 0.0;
299 
300  DBUG("PostNonLinearSolver crack volume computation.");
301 
302  const int process_id = 0;
303  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
304 
306  &LocalAssemblerInterface::computeCrackIntegral, _local_assemblers,
307  pv.getActiveElementIDs(), *_local_to_global_index_map, x,
308  _process_data.crack_volume);
309 
310  INFO("Integral of crack: {:g}", _process_data.crack_volume);
311 
313 }
314 
317 
318 } // namespace SmallDeformationNonlocal
319 } // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:34
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:29
Global vector based on Eigen vector.
Definition: EigenVector.h:28
void copyValues(std::vector< double > &u) const
Definition: EigenVector.h:106
bool isAxiallySymmetric() const
Definition: Mesh.h:131
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:76
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:103
Properties & getProperties()
Definition: Mesh.h:128
bool existsPropertyVector(std::string const &name) const
virtual void postTimestep(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt)
virtual void initialize(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
std::vector< std::size_t > const & getActiveElementIDs() const
std::vector< std::unique_ptr< IntegrationPointWriter > > _integration_point_writer
Definition: Process.h:354
SecondaryVariableCollection _secondary_variables
Definition: Process.h:335
Handles configuration of several secondary variables from the project file.
void addSecondaryVariable(std::string const &internal_name, SecondaryVariableFunctions &&fcts)
void assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
void assembleConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
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)
std::vector< std::unique_ptr< LocalAssemblerInterface > > _local_assemblers
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
NumLib::IterationResult postIterationConcreteProcess(GlobalVector const &x) override
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, int const process_id) override
void preAssemble(const std::size_t mesh_item_id, LocalAssemblerInterface &local_assembler, const NumLib::LocalToGlobalIndexMap &dof_table, const double t, double const dt, const GlobalVector &x)
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, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac)
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, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
IterationResult
Status flags telling the NonlinearSolver if an iteration succeeded.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:48
static const double t
@ BY_LOCATION
Ordering data by spatial location.
void createLocalAssemblers(std::vector< MeshLib::Element * > const &mesh_elements, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< std::unique_ptr< LocalAssemblerInterface >> &local_assemblers, NumLib::IntegrationOrder const integration_order, ExtraCtorArgs &&... extra_ctor_args)
IntegrationPointMetaData getIntegrationPointMetaData(MeshLib::Properties const &properties, std::string const &name)
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeMemberOnDereferenced(Method method, Container const &container, Args &&... args)