OGS
VectorMatrixAssembler.cpp
Go to the documentation of this file.
1
12
13#include <cassert>
14#include <functional> // for std::reference_wrapper.
15
20#include "Process.h"
21
22namespace 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 double const dt, 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, dt, 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, double const dt, std::vector<GlobalVector*> const& x,
46 std::vector<GlobalVector*> const& xdot, int const process_id,
48{
49 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
50 indices_of_processes.reserve(dof_tables.size());
51 transform(cbegin(dof_tables), cend(dof_tables),
52 back_inserter(indices_of_processes),
53 [&](auto const& dof_table)
54 { return NumLib::getIndices(mesh_item_id, dof_table); });
55
56 auto const& indices = indices_of_processes[process_id];
57 _local_M_data.clear();
58 _local_K_data.clear();
59 _local_b_data.clear();
60
61 std::size_t const number_of_processes = x.size();
62 // Monolithic scheme
63 if (number_of_processes == 1)
64 {
65 auto const local_x = x[process_id]->get(indices);
66 auto const local_xdot = xdot[process_id]->get(indices);
67 local_assembler.assemble(t, dt, local_x, local_xdot, _local_M_data,
69 }
70 else // Staggered scheme
71 {
72 auto local_coupled_xs =
73 getCoupledLocalSolutions(x, indices_of_processes);
74 auto const local_x = MathLib::toVector(local_coupled_xs);
75
76 auto local_coupled_xdots =
77 getCoupledLocalSolutions(xdot, indices_of_processes);
78 auto const local_xdot = MathLib::toVector(local_coupled_xdots);
79
80 local_assembler.assembleForStaggeredScheme(
81 t, dt, local_x, local_xdot, process_id, _local_M_data,
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, double const dt, std::vector<GlobalVector*> const& x,
111 std::vector<GlobalVector*> const& xdot, int const process_id,
113{
114 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
115 indices_of_processes.reserve(dof_tables.size());
116 transform(cbegin(dof_tables), cend(dof_tables),
117 back_inserter(indices_of_processes),
118 [&](auto const& dof_table)
119 { return NumLib::getIndices(mesh_item_id, dof_table); });
120
121 auto const& indices = indices_of_processes[process_id];
122
123 _local_M_data.clear();
124 _local_K_data.clear();
125 _local_b_data.clear();
126 _local_Jac_data.clear();
127
128 std::size_t const number_of_processes = x.size();
129 // Monolithic scheme
130 if (number_of_processes == 1)
131 {
132 auto const local_x = x[process_id]->get(indices);
133 auto const local_xdot = xdot[process_id]->get(indices);
134 _jacobian_assembler->assembleWithJacobian(
135 local_assembler, t, dt, local_x, local_xdot, _local_M_data,
137 }
138 else // Staggered scheme
139 {
140 auto local_coupled_xs =
141 getCoupledLocalSolutions(x, indices_of_processes);
142 auto const local_x = MathLib::toVector(local_coupled_xs);
143
144 auto local_coupled_xdots =
145 getCoupledLocalSolutions(xdot, indices_of_processes);
146 auto const local_xdot = MathLib::toVector(local_coupled_xdots);
147
148 _jacobian_assembler->assembleWithJacobianForStaggeredScheme(
149 local_assembler, t, dt, local_x, local_xdot, process_id,
151 }
152
153 auto const num_r_c = indices.size();
154 auto const r_c_indices =
156
157 if (!_local_M_data.empty())
158 {
159 auto const local_M = MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
160 M.add(r_c_indices, local_M);
161 }
162 if (!_local_K_data.empty())
163 {
164 auto const local_K = MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
165 K.add(r_c_indices, local_K);
166 }
167 if (!_local_b_data.empty())
168 {
169 assert(_local_b_data.size() == num_r_c);
170 b.add(indices, _local_b_data);
171 }
172 if (!_local_Jac_data.empty())
173 {
174 auto const local_Jac =
175 MathLib::toMatrix(_local_Jac_data, num_r_c, num_r_c);
176 Jac.add(r_c_indices, local_Jac);
177 }
178 else
179 {
180 OGS_FATAL(
181 "No Jacobian has been assembled! This might be due to programming "
182 "errors in the local assembler of the current process.");
183 }
184}
185
186} // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
int add(IndexType row, IndexType col, double val)
Definition: EigenMatrix.h:87
Global vector based on Eigen vector.
Definition: EigenVector.h:28
void add(IndexType rowId, double v)
add entry
Definition: EigenVector.h:79
double get(IndexType rowId) const
get entry
Definition: EigenVector.h:61
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
virtual void assembleForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
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::unique_ptr< AbstractJacobianAssembler > _jacobian_assembler
Used to assemble the Jacobian.
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)
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, 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 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)
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
static const double t
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:72
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
std::vector< double > getCoupledLocalSolutions(std::vector< GlobalVector * > const &global_solutions, std::vector< std::vector< GlobalIndexType > > const &indices)