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