v0.14.0
Schur.hpp
Go to the documentation of this file.
1 /**
2  * @file Schur.hpp
3  * @brief Assemble Schur complement
4  * @date 2023-02-02
5  *
6  * @copyright Copyright (c) 2023
7  *
8  * To create nested system of Schur complements, you push sequence of operator,
9  * to set up data on entities, and then assemble complements.
10  *
11  * \note Try separate floating points operations from book keeping. Also, align
12  * memory that blocks follow floating point operations.
13  *
14  */
15 
16 #ifndef __SCHUR_HPP__
17 #define __SCHUR_HPP__
18 
19 namespace MoFEM {
20 
21 constexpr const char MoFEM_BLOCK_MAT[] = "mofem_block_mat";
22 
23 /**
24  * @brief Structure to register events for Schur block assembly and solver
25  */
26 struct SchurEvents {
27  static PetscLogEvent MOFEM_EVENT_schurMatSetValues;
28  static PetscLogEvent MOFEM_EVENT_opSchurAssembleEnd;
30  static PetscLogEvent MOFEM_EVENT_BlockStructureMult;
31  static PetscLogEvent MOFEM_EVENT_BlockStructureSolve;
32  static PetscLogEvent MOFEM_EVENT_zeroRowsAndCols;
33  SchurEvents();
34 };
35 
36 struct SchurElemMats;
37 struct SchurElemMatsBlock;
39 
41 
42  OpSchurAssembleBase() = delete;
43 
44 private:
46 };
47 
49 
50 /**
51  * @brief Construct a new Op Schur Assemble End object
52  *
53  * @param fields_name list of fields (can be empty)
54  * @param field_ents list of entities on which schur complement is applied
55  * @param schur_aos maps dofs indices from main problem to schur complement
56  * @param schur_mat schur matrix
57  * @param sym_schur true if schur (matrix) complement is symmetric
58  * @param symm_op true if block diagonal is symmetric
59  */
61 createOpSchurAssembleEnd(std::vector<std::string> fields_name,
62  std::vector<boost::shared_ptr<Range>> field_ents,
65  bool sym_schur = false, bool symm_op = false);
66 
67 using SchurFieldPair = std::pair<std::string, std::string>;
68 
69 using SchurFEOpsFEandFields = std::vector<
70 
71  std::pair<std::string, std::vector<SchurFieldPair>>
72 
73  >;
74 
75 /**
76  * @brief Create a Mat Diag Blocks object
77  *
78  * @return Mat
79  */
80 boost::shared_ptr<BlockStructure> createBlockMatStructure(
81 
82  DM dm, //< dm
83  SchurFEOpsFEandFields schur_fe_op_vec //< block elements
84 
85 );
86 
87 using SchurShellMatData =
88  std::pair<SmartPetscObj<Mat>, boost::shared_ptr<BlockStructure>>;
89 
90 /**
91  * @brief Create a Schur Mat object
92  *
93  * @param dm
94  * @param data
95  * @return std::pair<SmartPetscObj<Mat>, boost::shared_ptr<BlockStructure>>
96  */
97 SchurShellMatData createBlockMat(DM dm, boost::shared_ptr<BlockStructure> data);
98 
99 using NestSchurData = std::tuple<
100 
101  std::array<SmartPetscObj<Mat>, 4>,
102  std::array<boost::shared_ptr<BlockStructure>, 4>,
103  boost::shared_ptr<BlockStructure>,
104  std::pair<SmartPetscObj<IS>, SmartPetscObj<IS>>
105 
106  >;
107 
108 /**
109  * @brief Get the Schur Nest Mat Array object
110  *
111  * @param dms schur dm, and block dm
112  * @param block mat A data
113  * @param fields_name list of fields
114  * @param field_ents list of entities on which schur complement is applied
115  * @param add_preconditioner_block add block for preconditioner
116  * @return boost::shared_ptr<NestSchurData>
117  */
118 boost::shared_ptr<NestSchurData> createSchurNestedMatrixStruture(
119  std::pair<SmartPetscObj<DM>, SmartPetscObj<DM>> dms,
120  boost::shared_ptr<BlockStructure> block_mat_data,
121 
122  std::vector<std::string> fields_name, //< a00 fields
123  std::vector<boost::shared_ptr<Range>> field_ents, //< a00 ents
124  bool add_preconditioner_block = false
125 
126 );
127 
128 /**
129  * @brief Switch preconditioner
130  *
131  * @param block_mat_data
132  * @return MoFEMErrorCode
133  */
135 schurSwitchPreconditioner(boost::shared_ptr<BlockStructure> block_mat_data);
136 
137 /**
138  * @brief Save block matrix as a mesh
139  *
140  * @param block_mat_data
141  * @param filename
142  * @return MoFEMErrorCode
143  */
145 schurSaveBlockMesh(boost::shared_ptr<BlockStructure> block_mat_data,
146  std::string filename);
147 
148 /**
149  * @brief Create a Mat Diag Blocks object
150  *
151  * \code {.cpp}
152  *
153  * auto [nested_mat, nested_data_ptr] = createSchurNestedMatrix(
154  *
155  * createSchurNestedMatrixStruture(
156  *
157  * {schur_dm, block_dm}, shell_data,
158  *
159  * {"TENSOR"}, {nullptr}
160  *
161  * )
162  *
163  * );
164  *
165  * \endcode
166  *
167  * @return Mat
168  */
169 std::pair<SmartPetscObj<Mat>, boost::shared_ptr<NestSchurData>>
170 createSchurNestedMatrix(boost::shared_ptr<NestSchurData> schur_net_data_ptr);
171 
172 template <>
173 MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr<NestSchurData>);
174 
175 /**
176  * @brief Set PC for A00 block
177  *
178  * @param pc
179  * @return MoFEMErrorCode
180  */
182 
184 
186  SchurBackendMatSetValuesPtr() = delete;
187  using MatSetValuesPtr = boost::function<MoFEMErrorCode(
188  Mat M, const EntitiesFieldData::EntData &row_data,
189  const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
190  InsertMode iora)>;
191  static MatSetValuesPtr matSetValuesPtr; ///< backend assembly function
192  static MatSetValuesPtr
193  matSetValuesBlockPtr; ///< backend assembly block mat function
194  static MatSetValuesPtr
195  matSetValuesPreconditionedBlockPtr; ///< backend assembly block
196  ///< preconditioner mat function
197 };
198 
199 template <>
202  const EntitiesFieldData::EntData &col_data,
203  const MatrixDouble &mat, InsertMode iora);
204 
205 template <>
206 inline MoFEMErrorCode
208  const VectorDouble &nf, InsertMode iora) {
209  return VecSetValues<EssentialBcStorage>(V, data, nf, iora);
210 }
211 
212 template <>
213 MoFEMErrorCode MatSetValues<AssemblyTypeSelector<SCHUR>>(
214  Mat M, const EntitiesFieldData::EntData &row_data,
215  const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
216  InsertMode iora);
217 
218 template <>
219 inline MoFEMErrorCode VecSetValues<AssemblyTypeSelector<SCHUR>>(
220  Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf,
221  InsertMode iora) {
222  return VecSetValues<SchurElemMats>(V, data, nf, iora);
223 }
224 
225 /***
226  * @brief Specialization of MatSetValues for BlockStructure
227  */
228 template <>
231  const EntitiesFieldData::EntData &col_data,
232  const MatrixDouble &mat, InsertMode iora);
233 
234 /***
235  * @brief Specialisation of MatSetValues for AssemblyTypeSelector<BLOCK_MAT>
236  */
237 template <>
238 inline MoFEMErrorCode MatSetValues<AssemblyTypeSelector<BLOCK_MAT>>(
239  Mat M, const EntitiesFieldData::EntData &row_data,
240  const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
241  InsertMode iora) {
242  return MatSetValues<BlockStructure>(M, row_data, col_data, mat, iora);
243 }
244 
245 /***
246  * @brief Specialisation of VecSetValues for AssemblyTypeSelector<BLOCK_MAT>
247  */
248 template <>
249 inline MoFEMErrorCode VecSetValues<AssemblyTypeSelector<BLOCK_MAT>>(
250  Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf,
251  InsertMode iora) {
252  return VecSetValues<EssentialBcStorage>(V, data, nf, iora);
253 }
254 
255 /***
256  * @brief Specialization of MatSetValues for SchurElemMatsBlock
257  */
258 template <>
261  const EntitiesFieldData::EntData &row_data,
262  const EntitiesFieldData::EntData &col_data,
263  const MatrixDouble &mat, InsertMode iora);
264 
265 /***
266  * @brief Specialization of VecSetValues for SchurElemMatsBlock
267  */
268 template <>
269 inline MoFEMErrorCode
271  const VectorDouble &nf, InsertMode iora) {
272  return VecSetValues<EssentialBcStorage>(V, data, nf, iora);
273 }
274 
275 /***
276  * @brief Specialisation of MatSetValues for AssemblyTypeSelector<BLOCK_SCHUR>
277  */
278 template <>
279 inline MoFEMErrorCode MatSetValues<AssemblyTypeSelector<BLOCK_SCHUR>>(
280  Mat M, const EntitiesFieldData::EntData &row_data,
281  const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
282  InsertMode iora) {
283  return MatSetValues<SchurElemMatsBlock>(M, row_data, col_data, mat, iora);
284 }
285 
286 /***
287  * @brief Specialisation of VecSetValues for AssemblyTypeSelector<BLOCK_SCHUR>
288  */
289 template <>
290 inline MoFEMErrorCode VecSetValues<AssemblyTypeSelector<BLOCK_SCHUR>>(
291  Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf,
292  InsertMode iora) {
293  return VecSetValues<EssentialBcStorage>(V, data, nf, iora);
294 }
295 
296 /***
297  * @brief Specialization of MatSetValues for SchurElemMatsPreconditionedBlock
298  */
299 template <>
301  Mat M, const EntitiesFieldData::EntData &row_data,
302  const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
303  InsertMode iora);
304 
305 /***
306  * @brief Specialisation of MatSetValues for
307  * AssemblyTypeSelector<BLOCK_PRECONDITIONER_SCHUR>
308  */
309 template <>
310 inline MoFEMErrorCode
311 MatSetValues<AssemblyTypeSelector<BLOCK_PRECONDITIONER_SCHUR>>(
312  Mat M, const EntitiesFieldData::EntData &row_data,
313  const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat,
314  InsertMode iora) {
315  return MatSetValues<SchurElemMatsPreconditionedBlock>(M, row_data, col_data,
316  mat, iora);
317 }
318 
319 /***
320  * @brief Specialisation of VecSetValues for AssemblyTypeSelector<BLOCK_SCHUR>
321  */
322 template <>
323 inline MoFEMErrorCode
324 VecSetValues<AssemblyTypeSelector<BLOCK_PRECONDITIONER_SCHUR>>(
325  Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf,
326  InsertMode iora) {
327  return ::VecSetValues(V, data.getIndices().size(),
328  &*data.getIndices().begin(), &*nf.begin(), iora);
329 }
330 
331 /**
332  * @deprecated do not use
333  */
334 DEPRECATED OpSchurAssembleBase *createOpSchurAssembleEnd(
335  std::vector<std::string> fields_name,
336  std::vector<boost::shared_ptr<Range>> field_ents,
337  std::vector<SmartPetscObj<AO>> sequence_of_aos,
338  std::vector<SmartPetscObj<Mat>> sequence_of_mats,
339  std::vector<bool> sym_schur, bool symm_op,
340  boost::shared_ptr<BlockStructure> diag_blocks = nullptr);
341 
342 /**
343  * @deprecated do not use
344  */
345 DEPRECATED OpSchurAssembleBase *createOpSchurAssembleEnd(
346  std::vector<std::string> fields_name,
347  std::vector<boost::shared_ptr<Range>> field_ents,
348  std::vector<SmartPetscObj<AO>> sequence_of_aos,
349  std::vector<SmartPetscObj<Mat>> sequence_of_mats,
350  std::vector<bool> sym_schur, std::vector<double> diag_eps, bool symm_op,
351  boost::shared_ptr<BlockStructure> diag_blocks = nullptr);
352 
353 /** @deprecated do not use, use createSchurNestedMatrixStruture instead */
354 DEPRECATED inline boost::shared_ptr<NestSchurData>
356  boost::shared_ptr<BlockStructure> block_mat_data,
357 
358  std::vector<std::string> fields_name, //< a00 fields
359  std::vector<boost::shared_ptr<Range>> field_ents, //< a00 ents
360  bool add_preconditioner_block = false
361 
362 ) {
363  return createSchurNestedMatrixStruture(dms, block_mat_data, fields_name,
364  field_ents, add_preconditioner_block);
365 }
366 
367 } // namespace MoFEM
368 
369 #endif //__SCHUR_HPP__
UBlasMatrix< double >
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
DEPRECATED
#define DEPRECATED
Definition: definitions.h:17
MoFEM::MatSetValues< BlockStructure >
MoFEMErrorCode MatSetValues< BlockStructure >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
Definition: Schur.cpp:2293
MoFEM::MoFEM_BLOCK_MAT
constexpr const char MoFEM_BLOCK_MAT[]
Definition: Schur.hpp:21
MoFEM::SchurFieldPair
std::pair< std::string, std::string > SchurFieldPair
Definition: Schur.hpp:67
MoFEM::OpSchurAssembleBase::OpSchurAssembleBase
OpSchurAssembleBase()=delete
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::SchurBackendMatSetValuesPtr::MatSetValuesPtr
boost::function< MoFEMErrorCode(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)> MatSetValuesPtr
Definition: Schur.hpp:190
MoFEM::SchurEvents::SchurEvents
SchurEvents()
Definition: Schur.cpp:299
MoFEM::SchurEvents::MOFEM_EVENT_BlockStructureMult
static PetscLogEvent MOFEM_EVENT_BlockStructureMult
Definition: Schur.hpp:30
MoFEM::SchurEvents::MOFEM_EVENT_schurMatSetValues
static PetscLogEvent MOFEM_EVENT_schurMatSetValues
Definition: Schur.hpp:27
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
MoFEM::createOpSchurAssembleEnd
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range >> field_ents, SmartPetscObj< AO > ao, SmartPetscObj< Mat > schur, bool sym_schur, bool symm_op)
Construct a new Op Schur Assemble End object.
Definition: Schur.cpp:2186
MoFEM::VecSetValues< SchurElemMatsBlock >
MoFEMErrorCode VecSetValues< SchurElemMatsBlock >(Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf, InsertMode iora)
Definition: Schur.hpp:270
MoFEM::createBlockMat
SchurShellMatData createBlockMat(DM dm, boost::shared_ptr< BlockStructure > data)
Create a Schur Mat object.
Definition: Schur.cpp:1379
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::SchurEvents::MOFEM_EVENT_BlockStructureSolve
static PetscLogEvent MOFEM_EVENT_BlockStructureSolve
Definition: Schur.hpp:31
MoFEM::schurSwitchPreconditioner
MoFEMErrorCode schurSwitchPreconditioner(boost::shared_ptr< BlockStructure > block_mat_data)
Switch preconditioner.
Definition: Schur.cpp:2365
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MoFEM::getNestSchurData
DEPRECATED boost::shared_ptr< NestSchurData > getNestSchurData(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM >> dms, boost::shared_ptr< BlockStructure > block_mat_data, std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range >> field_ents, bool add_preconditioner_block=false)
Definition: Schur.hpp:355
MoFEM::SchurBackendMatSetValuesPtr::SchurBackendMatSetValuesPtr
SchurBackendMatSetValuesPtr()=delete
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition: Schur.cpp:2181
MoFEM::SchurEvents::MOFEM_EVENT_opSchurAssembleEnd
static PetscLogEvent MOFEM_EVENT_opSchurAssembleEnd
Definition: Schur.hpp:28
MoFEM::SchurElemMatsBlock
Definition: Schur.cpp:2299
MoFEM::createSchurNestedMatrix
std::pair< SmartPetscObj< Mat >, boost::shared_ptr< NestSchurData > > createSchurNestedMatrix(boost::shared_ptr< NestSchurData > schur_net_data_ptr)
Create a Mat Diag Blocks object.
Definition: Schur.cpp:2153
MoFEM::MatSetValues< SchurElemMatsPreconditionedBlock >
MoFEMErrorCode MatSetValues< SchurElemMatsPreconditionedBlock >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
Definition: Schur.cpp:2356
MoFEM::SchurBackendMatSetValuesPtr::matSetValuesPtr
static MatSetValuesPtr matSetValuesPtr
backend assembly function
Definition: Schur.hpp:191
MoFEM::SchurElemMatsPreconditionedBlock
Definition: Schur.cpp:2332
MoFEM::SchurBackendMatSetValuesPtr
Definition: Schur.hpp:185
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
MoFEM::OpSchurAssembleBase
Definition: Schur.hpp:40
MoFEM::setSchurMatSolvePC
MoFEMErrorCode setSchurMatSolvePC(SmartPetscObj< PC > pc)
Definition: Schur.cpp:2241
MoFEM::SchurElemMats
Schur complement data storage.
Definition: Schur.cpp:147
MoFEM::setSchurA00MatSolvePC
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
Definition: Schur.cpp:2223
MoFEM::SchurBackendMatSetValuesPtr::matSetValuesPreconditionedBlockPtr
static MatSetValuesPtr matSetValuesPreconditionedBlockPtr
Definition: Schur.hpp:195
MoFEM::MatSetValues< SchurElemMats >
MoFEMErrorCode MatSetValues< SchurElemMats >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
Definition: Schur.cpp:2249
MoFEM::createSchurNestedMatrixStruture
boost::shared_ptr< NestSchurData > createSchurNestedMatrixStruture(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM >> dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range >> field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
Definition: Schur.cpp:1944
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::SchurEvents
Structure to register events for Schur block assembly and solver.
Definition: Schur.hpp:26
MoFEM::SchurShellMatData
std::pair< SmartPetscObj< Mat >, boost::shared_ptr< BlockStructure > > SchurShellMatData
Definition: Schur.hpp:88
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEM::DMMoFEMSetNestSchurData
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
Definition: DMMoFEM.cpp:1562
MoFEM::VecSetValues< SchurElemMats >
MoFEMErrorCode VecSetValues< SchurElemMats >(Vec V, const EntitiesFieldData::EntData &data, const VectorDouble &nf, InsertMode iora)
Definition: Schur.hpp:207
MoFEM::createBlockMatStructure
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
Definition: Schur.cpp:1009
MoFEM::SchurBackendMatSetValuesPtr::matSetValuesBlockPtr
static MatSetValuesPtr matSetValuesBlockPtr
backend assembly block mat function
Definition: Schur.hpp:193
MoFEM::SmartPetscObj< AO >
MoFEM::MatSetValues< SchurElemMatsBlock >
MoFEMErrorCode MatSetValues< SchurElemMatsBlock >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const MatrixDouble &mat, InsertMode iora)
Definition: Schur.cpp:2325
MoFEM::NestSchurData
std::tuple< std::array< SmartPetscObj< Mat >, 4 >, std::array< boost::shared_ptr< BlockStructure >, 4 >, boost::shared_ptr< BlockStructure >, std::pair< SmartPetscObj< IS >, SmartPetscObj< IS > > > NestSchurData
Definition: Schur.hpp:106
MoFEM::SchurEvents::MOFEM_EVENT_zeroRowsAndCols
static PetscLogEvent MOFEM_EVENT_zeroRowsAndCols
Definition: Schur.hpp:32
MoFEM::schurSaveBlockMesh
MoFEMErrorCode schurSaveBlockMesh(boost::shared_ptr< BlockStructure > block_mat_data, std::string filename)
Save block matrix as a mesh.
Definition: Schur.cpp:2380
MoFEM::SchurEvents::MOFEM_EVENT_BlockStructureSetValues
static PetscLogEvent MOFEM_EVENT_BlockStructureSetValues
Definition: Schur.hpp:29
MoFEM::SchurFEOpsFEandFields
std::vector< std::pair< std::string, std::vector< SchurFieldPair > > > SchurFEOpsFEandFields
Definition: Schur.hpp:73
MoFEM::VecSetValues< EssentialBcStorage >
MoFEMErrorCode VecSetValues< EssentialBcStorage >(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Set values to vector in operator.
Definition: FormsIntegrators.cpp:76
MoFEM::ForcesAndSourcesCore::UserDataOperator::UserDataOperator
UserDataOperator(const FieldSpace space, const char type=OPSPACE, const bool symm=true)
Definition: ForcesAndSourcesCore.cpp:1977