OGS 6.2.0-405-gb717f6088
VectorMatrixAssembler.cpp
Go to the documentation of this file.
1 
10 #include "VectorMatrixAssembler.h"
11 
12 #include <cassert>
13 #include <functional> // for std::reference_wrapper.
14 
18 
20 #include "Process.h"
21 
22 namespace ProcessLib
23 {
25  std::unique_ptr<AbstractJacobianAssembler>&& jacobian_assembler)
26  : _jacobian_assembler(std::move(jacobian_assembler))
27 {
28 }
29 
31  const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
32  const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
33  const GlobalVector& x)
34 {
35  auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
36  auto const local_x = x.get(indices);
37 
38  local_assembler.preAssemble(t, local_x);
39 }
40 
42  const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
43  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
44  dof_tables,
45  const double t, const GlobalVector& x, GlobalMatrix& M, GlobalMatrix& K,
46  GlobalVector& b, CoupledSolutionsForStaggeredScheme const* const cpl_xs)
47 {
48  std::vector<std::vector<GlobalIndexType>> indices_of_processes;
49  indices_of_processes.reserve(dof_tables.size());
50  for (auto dof_table : dof_tables)
51  {
52  indices_of_processes.emplace_back(
53  NumLib::getIndices(mesh_item_id, dof_table.get()));
54  }
55 
56  auto const& indices = (cpl_xs == nullptr)
57  ? indices_of_processes[0]
58  : indices_of_processes[cpl_xs->process_id];
59  _local_M_data.clear();
60  _local_K_data.clear();
61  _local_b_data.clear();
62 
63  if (cpl_xs == nullptr)
64  {
65  auto const local_x = x.get(indices);
66  local_assembler.assemble(t, local_x, _local_M_data, _local_K_data,
68  }
69  else
70  {
71  auto local_coupled_xs0 =
72  getPreviousLocalSolutions(*cpl_xs, indices_of_processes);
73  auto local_coupled_xs =
74  getCurrentLocalSolutions(*cpl_xs, indices_of_processes);
75 
76  ProcessLib::LocalCoupledSolutions local_coupled_solutions(
77  cpl_xs->dt, cpl_xs->process_id, std::move(local_coupled_xs0),
78  std::move(local_coupled_xs));
79 
80  local_assembler.assembleForStaggeredScheme(t, _local_M_data,
82  local_coupled_solutions);
83  }
84 
85  auto const num_r_c = indices.size();
86  auto const r_c_indices =
88 
89  if (!_local_M_data.empty())
90  {
91  auto const local_M = MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
92  M.add(r_c_indices, local_M);
93  }
94  if (!_local_K_data.empty())
95  {
96  auto const local_K = MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
97  K.add(r_c_indices, local_K);
98  }
99  if (!_local_b_data.empty())
100  {
101  assert(_local_b_data.size() == num_r_c);
102  b.add(indices, _local_b_data);
103  }
104 }
105 
107  std::size_t const mesh_item_id, LocalAssemblerInterface& local_assembler,
108  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
109  dof_tables,
110  const double t, GlobalVector const& x, GlobalVector const& xdot,
111  const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
112  GlobalVector& b, GlobalMatrix& Jac,
113  CoupledSolutionsForStaggeredScheme const* const cpl_xs)
114 {
115  std::vector<std::vector<GlobalIndexType>> indices_of_processes;
116  indices_of_processes.reserve(dof_tables.size());
117  for (auto dof_table : dof_tables)
118  {
119  indices_of_processes.emplace_back(
120  NumLib::getIndices(mesh_item_id, dof_table.get()));
121  }
122 
123  auto const& indices = (cpl_xs == nullptr)
124  ? indices_of_processes[0]
125  : indices_of_processes[cpl_xs->process_id];
126  auto const local_xdot = xdot.get(indices);
127 
128  _local_M_data.clear();
129  _local_K_data.clear();
130  _local_b_data.clear();
131  _local_Jac_data.clear();
132 
133  if (cpl_xs == nullptr)
134  {
135  auto const local_x = x.get(indices);
136  _jacobian_assembler->assembleWithJacobian(
137  local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx,
139  }
140  else
141  {
142  auto local_coupled_xs0 =
143  getPreviousLocalSolutions(*cpl_xs, indices_of_processes);
144  auto local_coupled_xs =
145  getCurrentLocalSolutions(*cpl_xs, indices_of_processes);
146 
147  ProcessLib::LocalCoupledSolutions local_coupled_solutions(
148  cpl_xs->dt, cpl_xs->process_id, std::move(local_coupled_xs0),
149  std::move(local_coupled_xs));
150 
151  _jacobian_assembler->assembleWithJacobianForStaggeredScheme(
152  local_assembler, t, local_xdot, dxdot_dx, dx_dx, _local_M_data,
154  local_coupled_solutions);
155  }
156 
157  auto const num_r_c = indices.size();
158  auto const r_c_indices =
160 
161  if (!_local_M_data.empty())
162  {
163  auto const local_M = MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
164  M.add(r_c_indices, local_M);
165  }
166  if (!_local_K_data.empty())
167  {
168  auto const local_K = MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
169  K.add(r_c_indices, local_K);
170  }
171  if (!_local_b_data.empty())
172  {
173  assert(_local_b_data.size() == num_r_c);
174  b.add(indices, _local_b_data);
175  }
176  if (!_local_Jac_data.empty())
177  {
178  auto const local_Jac =
179  MathLib::toMatrix(_local_Jac_data, num_r_c, num_r_c);
180  Jac.add(r_c_indices, local_Jac);
181  }
182  else
183  {
184  OGS_FATAL(
185  "No Jacobian has been assembled! This might be due to programming "
186  "errors in the local assembler of the current process.");
187  }
188 }
189 
190 } // namespace ProcessLib
virtual void preAssemble(double const, std::vector< double > const &)
virtual void assembleForStaggeredScheme(double const t, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, LocalCoupledSolutions const &coupled_solutions)
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
void preAssemble(const std::size_t mesh_item_id, LocalAssemblerInterface &local_assembler, const NumLib::LocalToGlobalIndexMap &dof_table, const double t, const GlobalVector &x)
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)
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
VectorMatrixAssembler(std::unique_ptr< AbstractJacobianAssembler > &&jacobian_assembler)
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)
std::vector< std::vector< double > > getPreviousLocalSolutions(const CoupledSolutionsForStaggeredScheme &cpl_xs, const std::vector< std::vector< GlobalIndexType >> &indices)
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:71
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
virtual void assemble(double const t, std::vector< double > const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
std::vector< std::vector< double > > getCurrentLocalSolutions(const CoupledSolutionsForStaggeredScheme &cpl_xs, const std::vector< std::vector< GlobalIndexType >> &indices)
std::unique_ptr< AbstractJacobianAssembler > _jacobian_assembler
Used to assemble the Jacobian.