OGS
VectorMatrixAssembler.cpp
Go to the documentation of this file.
1
12
13#include <cassert>
14
19
20namespace ProcessLib
21{
23 AbstractJacobianAssembler& jacobian_assembler)
24 : _jacobian_assembler(jacobian_assembler)
25{
26}
27
29 const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
30 const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
31 double const dt, const GlobalVector& x)
32{
33 auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
34 auto const local_x = x.get(indices);
35
36 local_assembler.preAssemble(t, dt, local_x);
37}
38
40 const std::size_t mesh_item_id, LocalAssemblerInterface& local_assembler,
41 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
42 const double t, double const dt, std::vector<GlobalVector*> const& x,
43 std::vector<GlobalVector*> const& x_prev, int const process_id,
45{
46 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
47 indices_of_processes.reserve(dof_tables.size());
48 transform(cbegin(dof_tables), cend(dof_tables),
49 back_inserter(indices_of_processes),
50 [&](auto const dof_table)
51 { return NumLib::getIndices(mesh_item_id, *dof_table); });
52
53 auto const& indices = indices_of_processes[process_id];
54 _local_M_data.clear();
55 _local_K_data.clear();
56 _local_b_data.clear();
57
58 std::size_t const number_of_processes = x.size();
59 // Monolithic scheme
60 if (number_of_processes == 1)
61 {
62 auto const local_x = x[process_id]->get(indices);
63 auto const local_x_prev = x_prev[process_id]->get(indices);
64 local_assembler.assemble(t, dt, local_x, local_x_prev, _local_M_data,
66 }
67 else // Staggered scheme
68 {
69 auto local_coupled_xs =
70 getCoupledLocalSolutions(x, indices_of_processes);
71 auto const local_x = MathLib::toVector(local_coupled_xs);
72
73 auto local_coupled_x_prevs =
74 getCoupledLocalSolutions(x_prev, indices_of_processes);
75 auto const local_x_prev = MathLib::toVector(local_coupled_x_prevs);
76
77 local_assembler.assembleForStaggeredScheme(
78 t, dt, local_x, local_x_prev, process_id, _local_M_data,
80 }
81
82 auto const num_r_c = indices.size();
83 auto const r_c_indices =
85
86 if (M && !_local_M_data.empty())
87 {
88 auto const local_M = MathLib::toMatrix(_local_M_data, num_r_c, num_r_c);
89 M->add(r_c_indices, local_M);
90 }
91 if (K && !_local_K_data.empty())
92 {
93 auto const local_K = MathLib::toMatrix(_local_K_data, num_r_c, num_r_c);
94 K->add(r_c_indices, local_K);
95 }
96 if (b && !_local_b_data.empty())
97 {
98 assert(_local_b_data.size() == num_r_c);
99 b->add(indices, _local_b_data);
100 }
101
102 _local_output(t, process_id, mesh_item_id, _local_M_data, _local_K_data,
104}
105
107 std::size_t const mesh_item_id, LocalAssemblerInterface& local_assembler,
108 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
109 const double t, double const dt, std::vector<GlobalVector*> const& x,
110 std::vector<GlobalVector*> const& x_prev, int const process_id,
111 GlobalVector* b, GlobalMatrix* Jac)
112{
113 std::vector<std::vector<GlobalIndexType>> indices_of_processes;
114 indices_of_processes.reserve(dof_tables.size());
115 transform(cbegin(dof_tables), cend(dof_tables),
116 back_inserter(indices_of_processes),
117 [&](auto const dof_table)
118 { return NumLib::getIndices(mesh_item_id, *dof_table); });
119
120 auto const& indices = indices_of_processes[process_id];
121
122 _local_b_data.clear();
123 _local_Jac_data.clear();
124
125 std::size_t const number_of_processes = x.size();
126 // Monolithic scheme
127 if (number_of_processes == 1)
128 {
129 auto const local_x = x[process_id]->get(indices);
130 auto const local_x_prev = x_prev[process_id]->get(indices);
132 local_assembler, t, dt, local_x, local_x_prev, _local_b_data,
134 }
135 else // Staggered scheme
136 {
137 auto local_coupled_xs =
138 getCoupledLocalSolutions(x, indices_of_processes);
139 auto const local_x = MathLib::toVector(local_coupled_xs);
140
141 auto local_coupled_x_prevs =
142 getCoupledLocalSolutions(x_prev, indices_of_processes);
143 auto const local_x_prev = MathLib::toVector(local_coupled_x_prevs);
144
146 local_assembler, t, dt, local_x, local_x_prev, process_id,
148 }
149
150 auto const num_r_c = indices.size();
151 auto const r_c_indices =
153
154 if (b && !_local_b_data.empty())
155 {
156 assert(_local_b_data.size() == num_r_c);
157 b->add(indices, _local_b_data);
158 }
159 if (Jac && !_local_Jac_data.empty())
160 {
161 auto const local_Jac =
162 MathLib::toMatrix(_local_Jac_data, num_r_c, num_r_c);
163 Jac->add(r_c_indices, local_Jac);
164 }
165 else
166 {
167 OGS_FATAL(
168 "No Jacobian has been assembled! This might be due to programming "
169 "errors in the local assembler of the current process.");
170 }
171
172 _local_output(t, process_id, mesh_item_id, _local_b_data, _local_Jac_data);
173}
174
175} // 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:25
void add(IndexType rowId, double v)
add entry
Definition EigenVector.h:76
double get(IndexType rowId) const
get entry
Definition EigenVector.h:58
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
Base class for Jacobian assemblers.
virtual void assembleWithJacobian(LocalAssemblerInterface &local_assembler, double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)=0
virtual void assembleWithJacobianForStaggeredScheme(LocalAssemblerInterface &, double const, double const, Eigen::VectorXd const &, Eigen::VectorXd const &, int const, std::vector< double > &, std::vector< double > &)
virtual void preAssemble(double const, double const, std::vector< double > const &)
virtual void assembleForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
virtual void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix *M, GlobalMatrix *K, GlobalVector *b)
VectorMatrixAssembler(AbstractJacobianAssembler &jacobian_assembler)
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)
Assembly::LocalMatrixOutput _local_output
AbstractJacobianAssembler & _jacobian_assembler
Used to assemble the Jacobian.
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalVector *b, GlobalMatrix *Jac)
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
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)