v0.14.0
schur_test_diag_mat.cpp
Go to the documentation of this file.
1 /**
2  * @file schur_test_diag_mat.cpp
3  * @example schur_test_diag_mat.cpp
4  * @brief Test block matrix and Schur complement matrix.
5  *
6  * @copyright Copyright (c) 2024
7  *
8  */
9 
10 #include <MoFEM.hpp>
11 using namespace MoFEM;
12 
13 static char help[] = "...\n\n";
14 
15 template <int DIM> struct ElementsAndOps {};
16 
17 template <> struct ElementsAndOps<2> {
18  using DomainEle = PipelineManager::FaceEle;
19 };
20 
21 template <> struct ElementsAndOps<3> {
22  using DomainEle = VolumeElementForcesAndSourcesCore;
23 };
24 
25 //! [Define dimension]
26 constexpr int SPACE_DIM = 3;
27 constexpr int FIELD_DIM = SPACE_DIM;
28 
29 constexpr AssemblyType A =
30  AssemblyType::BLOCK_MAT; //< selected assembly type
31 constexpr IntegrationType I =
32  IntegrationType::GAUSS; //< selected integration type
33 
34 using EntData = EntitiesFieldData::EntData;
35 using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
36 using DomainEleOp = DomainEle::UserDataOperator;
37 
38 SmartPetscObj<Mat> petsc_mat;
39 SmartPetscObj<Mat> block_mat;
40 
41 template <>
42 typename MoFEM::OpBaseImpl<PETSC, DomainEleOp>::MatSetValuesHook
43  MoFEM::OpBaseImpl<PETSC, DomainEleOp>::matSetValuesHook =
44  [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
45  const EntitiesFieldData::EntData &row_data,
46  const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
47  return MatSetValues<AssemblyTypeSelector<PETSC>>(
48  petsc_mat, row_data, col_data, m, ADD_VALUES);
49  };
50 
51 template <>
52 typename MoFEM::OpBaseImpl<BLOCK_MAT, DomainEleOp>::MatSetValuesHook
53  MoFEM::OpBaseImpl<BLOCK_MAT, DomainEleOp>::matSetValuesHook =
54  [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
55  const EntitiesFieldData::EntData &row_data,
56  const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
57  return MatSetValues<AssemblyTypeSelector<BLOCK_MAT>>(
58  block_mat, row_data, col_data, m, ADD_VALUES);
59  };
60 
61 template <>
62 typename MoFEM::OpBaseImpl<BLOCK_SCHUR, DomainEleOp>::MatSetValuesHook
63  MoFEM::OpBaseImpl<BLOCK_SCHUR, DomainEleOp>::matSetValuesHook =
64  [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
65  const EntitiesFieldData::EntData &row_data,
66  const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
67  return MatSetValues<AssemblyTypeSelector<BLOCK_SCHUR>>(
68  block_mat, row_data, col_data, m, ADD_VALUES);
69  };
70 
71 template <>
72 typename MoFEM::OpBaseImpl<BLOCK_PRECONDITIONER_SCHUR,
73  DomainEleOp>::MatSetValuesHook
74  MoFEM::OpBaseImpl<BLOCK_PRECONDITIONER_SCHUR,
75  DomainEleOp>::matSetValuesHook =
76  [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
77  const EntitiesFieldData::EntData &row_data,
78  const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
79  return MatSetValues<AssemblyTypeSelector<BLOCK_PRECONDITIONER_SCHUR>>(
80  block_mat, row_data, col_data, m, ADD_VALUES);
81  };
82 
83 constexpr bool debug = false;
84 
85 int main(int argc, char *argv[]) {
86 
87  // initialize petsc
88  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
89 
90  try {
91 
92  // Create MoAB database
93  moab::Core moab_core;
94  moab::Interface &moab = moab_core;
95 
96  // Create MoFEM database and link it to MoAB
97  MoFEM::Core mofem_core(moab);
98  MoFEM::Interface &m_field = mofem_core;
99 
100  // Register DM Manager
101  DMType dm_name = "DMMOFEM";
102  CHKERR DMRegister_MoFEM(dm_name);
103 
104  // Add logging channel for example
105  auto core_log = logging::core::get();
106  core_log->add_sink(
107  LogManager::createSink(LogManager::getStrmWorld(), "Timeline"));
108  LogManager::setLog("Timeline");
109  MOFEM_LOG_TAG("Timeline", "Timeline");
110 
111  // Simple interface
112  auto simple = m_field.getInterface<Simple>();
113 
114  // get options from command line
115  CHKERR simple->getOptions();
116  // load mesh file
117  CHKERR simple->loadFile();
118 
119  CHKERR simple->addDomainField("V", H1, AINSWORTH_LEGENDRE_BASE, FIELD_DIM);
120  CHKERR simple->addDomainField("T", L2, AINSWORTH_LEGENDRE_BASE, FIELD_DIM);
121  CHKERR simple->addDomainField("S", L2, AINSWORTH_LEGENDRE_BASE, FIELD_DIM);
122  CHKERR simple->addDomainField("O", L2, AINSWORTH_LEGENDRE_BASE, FIELD_DIM);
123 
124  // set fields order, i.e. for most first cases order is sufficient.
125  int order = 3;
126  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
127  CHKERR simple->setFieldOrder("V", order);
128  CHKERR simple->setFieldOrder("T", order);
129  CHKERR simple->setFieldOrder("S", order - 1);
130  CHKERR simple->setFieldOrder("O", order - 2);
131 
132  // setup problem
133  CHKERR simple->setUp();
134 
135  // If block "VOL" exist, first index is removed from "T" field. That test
136  // blocks handling when some DOFs are removed.
137  auto bc_mng = m_field.getInterface<BcManager>();
138  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "VOL",
139  "T", 0, 1);
140 
141  auto schur_dm = createDM(m_field.get_comm(), "DMMOFEM");
142  CHKERR DMMoFEMCreateSubDM(schur_dm, simple->getDM(), "SCHUR");
143  CHKERR DMMoFEMSetSquareProblem(schur_dm, PETSC_TRUE);
144  CHKERR DMMoFEMAddElement(schur_dm, simple->getDomainFEName());
145  CHKERR DMMoFEMAddSubFieldRow(schur_dm, "V");
146  CHKERR DMMoFEMAddSubFieldCol(schur_dm, "V");
147  CHKERR DMSetUp(schur_dm);
148 
149  auto block_dm = createDM(m_field.get_comm(), "DMMOFEM");
150  CHKERR DMMoFEMCreateSubDM(block_dm, simple->getDM(), "BLOCK");
151  CHKERR DMMoFEMSetSquareProblem(block_dm, PETSC_TRUE);
152  CHKERR DMMoFEMAddElement(block_dm, simple->getDomainFEName());
153  CHKERR DMMoFEMAddSubFieldRow(block_dm, "T");
154  CHKERR DMMoFEMAddSubFieldCol(block_dm, "T");
155  CHKERR DMMoFEMAddSubFieldRow(block_dm, "S");
156  CHKERR DMMoFEMAddSubFieldCol(block_dm, "S");
157  CHKERR DMMoFEMAddSubFieldRow(block_dm, "O");
158  CHKERR DMMoFEMAddSubFieldCol(block_dm, "O");
159  CHKERR DMSetUp(block_dm);
160 
161  petsc_mat = createDMMatrix(simple->getDM());
162  auto S = createDMMatrix(schur_dm);
163 
164  auto shell_data =
165  createBlockMat(simple->getDM(),
166 
167  createBlockMatStructure(simple->getDM(),
168 
169  {{simple->getDomainFEName(),
170 
171  {{"V", "V"},
172  {"T", "T"},
173  {"S", "S"},
174  {"O", "O"},
175  {"V", "T"},
176  {"T", "V"},
177  {"S", "T"},
178  {"T", "S"},
179  {"T", "O"},
180  {"O", "T"},
181  {"S", "O"},
182  {"O", "S"}
183 
184  }}}
185 
186  )
187 
188  );
189 
190  auto [mat, block_data_ptr] = shell_data;
191  block_mat = mat;
192 
193  std::vector<std::string> fields{"T", "S", "O"};
194 
195  auto [nested_mat, nested_data_ptr] = createSchurNestedMatrix(
196 
197  createSchurNestedMatrixStruture(
198 
199  {schur_dm, block_dm}, block_data_ptr,
200 
201  fields, {nullptr, nullptr, nullptr}, true
202 
203  )
204 
205  );
206 
207  using OpMassPETSCAssemble = FormsIntegrators<DomainEleOp>::Assembly<
208  PETSC>::BiLinearForm<GAUSS>::OpMass<1, FIELD_DIM>;
209  using OpMassBlockAssemble = FormsIntegrators<DomainEleOp>::Assembly<
210  BLOCK_SCHUR>::BiLinearForm<GAUSS>::OpMass<1, FIELD_DIM>;
211  using OpMassBlockPreconditionerAssemble =
212  FormsIntegrators<DomainEleOp>::Assembly<BLOCK_PRECONDITIONER_SCHUR>::
213  BiLinearForm<GAUSS>::OpMass<1, FIELD_DIM>;
214 
215  // get operators tester
216  auto pip_mng = m_field.getInterface<PipelineManager>(); // get interface to
217  // pipeline manager
218 
219  auto close_zero = [](double, double, double) { return 1; };
220  auto beta = [](double, double, double) { return -1./2; };
221  auto gamma = [](double, double, double) { return -1. / 4; };
222 
223  pip_mng->setDomainLhsIntegrationRule(
224  [](int, int, int o) { return 2 * o; });
225  auto &pip_lhs = pip_mng->getOpDomainLhsPipeline();
226 
227  pip_lhs.push_back(new OpMassPETSCAssemble("V", "V"));
228  // pip_lhs.push_back(new OpMassPETSCAssemble("T", "T"));
229  pip_lhs.push_back(new OpMassPETSCAssemble("V", "T"));
230  pip_lhs.push_back(new OpMassPETSCAssemble("T", "V"));
231  pip_lhs.push_back(new OpMassPETSCAssemble("S", "S", close_zero));
232  pip_lhs.push_back(new OpMassPETSCAssemble("S", "T", beta));
233  pip_lhs.push_back(new OpMassPETSCAssemble("T", "S", beta));
234  pip_lhs.push_back(new OpMassPETSCAssemble("O", "O", close_zero));
235  pip_lhs.push_back(new OpMassPETSCAssemble("T", "O", beta));
236  pip_lhs.push_back(new OpMassPETSCAssemble("O", "T", beta));
237  pip_lhs.push_back(new OpMassPETSCAssemble("S", "O", gamma));
238  pip_lhs.push_back(new OpMassPETSCAssemble("O", "S", gamma));
239 
240  pip_lhs.push_back(createOpSchurAssembleBegin());
241  pip_lhs.push_back(new OpMassBlockAssemble("V", "V"));
242  // pip_lhs.push_back(new OpMassBlockAssemble("T", "T"));
243  pip_lhs.push_back(new OpMassBlockAssemble("V", "T"));
244  pip_lhs.push_back(new OpMassBlockAssemble("T", "V"));
245  pip_lhs.push_back(new OpMassBlockAssemble("S", "S", close_zero));
246  pip_lhs.push_back(new OpMassBlockAssemble("S", "T", beta));
247  pip_lhs.push_back(new OpMassBlockAssemble("T", "S", beta));
248  pip_lhs.push_back(new OpMassBlockAssemble("O", "O", close_zero));
249  pip_lhs.push_back(new OpMassBlockAssemble("T", "O", beta));
250  pip_lhs.push_back(new OpMassBlockAssemble("O", "T", beta));
251  pip_lhs.push_back(new OpMassBlockAssemble("S", "O", gamma));
252  pip_lhs.push_back(new OpMassBlockAssemble("O", "S", gamma));
253  pip_lhs.push_back(new OpMassBlockPreconditionerAssemble("T", "T"));
254 
255  auto schur_is = getDMSubData(schur_dm)->getSmartRowIs();
256  auto ao_up = createAOMappingIS(schur_is, PETSC_NULL);
257 
258  pip_lhs.push_back(createOpSchurAssembleEnd(
259 
260  fields, {nullptr, nullptr, nullptr}, ao_up, S, true, true)
261 
262  );
263 
264  {
265  MOFEM_LOG_CHANNEL("Timeline");
266  MOFEM_LOG_TAG("Timeline", "timer");
267  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
268  MOFEM_LOG("Timeline", Sev::inform) << "Assemble start";
269  CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
270  simple->getDomainFEName(),
271  pip_mng->getDomainLhsFE());
272  MOFEM_LOG("Timeline", Sev::inform) << "Assemble end";
273  }
274  {
275  MOFEM_LOG_CHANNEL("Timeline");
276  MOFEM_LOG_TAG("Timeline", "timer");
277  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
278  MOFEM_LOG("Timeline", Sev::inform) << "Mat assemble start";
279  CHKERR MatAssemblyBegin(petsc_mat, MAT_FINAL_ASSEMBLY);
280  CHKERR MatAssemblyEnd(petsc_mat, MAT_FINAL_ASSEMBLY);
281  MOFEM_LOG("Timeline", Sev::inform) << "Mat assemble end";
282  }
283 
284  auto get_random_vector = [&](auto dm) {
285  auto v = createDMVector(dm);
286  PetscRandom rctx;
287  PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
288  CHK_MOAB_THROW(VecSetRandom(v, rctx), "generate rand vector");
289  PetscRandomDestroy(&rctx);
290  return v;
291  };
292  auto v = get_random_vector(simple->getDM());
293 
294  auto y_petsc = createDMVector(simple->getDM());
295  auto y_block = createDMVector(simple->getDM());
296 
297  auto test = [](auto msg, auto y, double norm0) {
298  MoFEMFunctionBegin;
299 
300  double eps = 1e-10;
301  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-eps", &eps, PETSC_NULL);
302 
303  PetscReal norm;
304  CHKERR VecNorm(y, NORM_2, &norm);
305  norm = norm / norm0;
306  MOFEM_LOG_CHANNEL("WORLD");
307  MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "TestBlockMat")
308  << msg << ": norm of difference: " << norm;
309  if (norm > eps || std::isnan(norm) || std::isinf(norm)) {
310  SETERRQ(PETSC_COMM_WORLD, 1, "norm of difference is too big");
311  }
312  MoFEMFunctionReturn(0);
313  };
314 
315  std::vector<int> zero_rows_and_cols = {
316  0, 1, 10, 20,
317  500}; // not to remove dofs for TENSOR filed, inverse will not work
318  CHKERR MatZeroRowsColumns(petsc_mat, zero_rows_and_cols.size(),
319  &*zero_rows_and_cols.begin(), 1, PETSC_NULL,
320  PETSC_NULL);
321  CHKERR MatZeroRowsColumns(block_mat, zero_rows_and_cols.size(),
322  &*zero_rows_and_cols.begin(), 1, PETSC_NULL,
323  PETSC_NULL);
324 
325  {
326  MOFEM_LOG_CHANNEL("Timeline");
327  MOFEM_LOG_TAG("Timeline", "timer");
328  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
329  MOFEM_LOG("Timeline", Sev::inform)
330  << "MatMult(petsc_mat, v, y_petsc) star";
331  CHKERR MatMult(petsc_mat, v, y_petsc);
332  MOFEM_LOG("Timeline", Sev::inform)
333  << "MatMult(petsc_mat, v, y_petsc) end";
334  }
335  double nrm0;
336  CHKERR VecNorm(y_petsc, NORM_2, &nrm0);
337 
338  {
339  MOFEM_LOG_CHANNEL("Timeline");
340  MOFEM_LOG_TAG("Timeline", "timer");
341  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
342  MOFEM_LOG("Timeline", Sev::inform)
343  << "MatMult(block_mat, v, y_block) star";
344  CHKERR MatMult(block_mat, v, y_block);
345  MOFEM_LOG("Timeline", Sev::inform)
346  << "MatMult(block_mat, v, y_block) end";
347  }
348 
349  CHKERR VecAXPY(y_petsc, -1.0, y_block);
350  CHKERR test("mult", y_petsc, nrm0);
351 
352  MOFEM_LOG_CHANNEL("Timeline");
353  MOFEM_LOG_TAG("Timeline", "timer");
354 
355  CHKERR MatMult(petsc_mat, v, y_petsc);
356  CHKERR MatMult(block_mat, v, y_block);
357  CHKERR MatMultAdd(petsc_mat, v, y_petsc, y_petsc);
358  CHKERR MatMultAdd(block_mat, v, y_block, y_block);
359  CHKERR VecAXPY(y_petsc, -1.0, y_block);
360 
361  CHKERR test("mult add", y_petsc, nrm0);
362  CHKERR schurSwitchPreconditioner(std::get<1>(*nested_data_ptr)[3]);
363  auto y_nested = createDMVector(simple->getDM());
364  CHKERR MatMult(petsc_mat, v, y_petsc);
365 
366  CHKERR MatMult(nested_mat, v, y_nested);
367  CHKERR VecAXPY(y_petsc, -1.0, y_nested);
368 
369  CHKERR test("mult nested", y_petsc, nrm0);
370 
371  CHKERR schurSwitchPreconditioner(std::get<1>(*nested_data_ptr)[3]);
372  auto diag_mat = std::get<0>(*nested_data_ptr)[3];
373  auto diag_block_x = get_random_vector(block_dm);
374  auto diag_block_f = createDMVector(block_dm);
375  auto block_solved_x = createDMVector(block_dm);
376  CHKERR MatMult(diag_mat, diag_block_x, diag_block_f);
377  // That is if one like to use MatSolve directly, not though PC, as it is
378  // below
379  // CHKERR MatSolve(diag_mat, diag_block_f, block_solved_x);
380 
381  // set matrix type to shell, set data
382  CHKERR DMSetMatType(block_dm, MATSHELL);
383  CHKERR DMMoFEMSetBlocMatData(block_dm, std::get<1>(*nested_data_ptr)[3]);
384  // set empty operator, since block data are already calculated
385  CHKERR DMKSPSetComputeOperators(
386  block_dm,
387  [](KSP, Mat, Mat, void *) {
388  MOFEM_LOG("WORLD", Sev::inform) << "empty operator";
389  return 0;
390  },
391  nullptr);
392 
393  auto ksp = createKSP(m_field.get_comm());
394  CHKERR KSPSetDM(ksp, block_dm);
395  CHKERR KSPSetFromOptions(ksp);
396 
397  // set preconditioner to block mat
398  auto get_pc = [](auto ksp) {
399  PC pc_raw;
400  CHKERR KSPGetPC(ksp, &pc_raw);
401  return SmartPetscObj<PC>(pc_raw, true); // bump reference
402  };
403  CHKERR setSchurA00MatSolvePC(get_pc(ksp));
404  CHKERR KSPSetUp(ksp);
405 
406  CHKERR VecZeroEntries(block_solved_x);
407  CHKERR KSPSolve(ksp, diag_block_f, block_solved_x);
408 
409  auto diag_block_f_test = createDMVector(block_dm);
410  CHKERR MatMult(diag_mat, block_solved_x, diag_block_f_test);
411  CHKERR VecAXPY(diag_block_f_test, -1.0, diag_block_f);
412  CHKERR test("diag solve", diag_block_f_test, nrm0);
413 
414  if (m_field.get_comm_rank() == 0) {
415  CHKERR schurSaveBlockMesh(block_data_ptr, "block_mesh.vtk");
416  }
417 
418  petsc_mat.reset();
419  block_mat.reset();
420  }
421  CATCH_ERRORS;
422 
423  // finish work cleaning memory, getting statistics, etc.
424  MoFEM::Core::Finalize();
425 
426  return 0;
427 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
MoFEM::createAOMappingIS
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
Definition: PetscSmartObj.hpp:318
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
help
static char help[]
Definition: schur_test_diag_mat.cpp:13
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:280
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
MoFEM::BLOCK_PRECONDITIONER_SCHUR
@ BLOCK_PRECONDITIONER_SCHUR
Definition: FormsIntegrators.hpp:109
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::BLOCK_MAT
@ BLOCK_MAT
Definition: FormsIntegrators.hpp:107
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:261
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::getDMSubData
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition: DMMoFEM.hpp:1157
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:1056
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
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
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
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::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
MoFEM::schurSwitchPreconditioner
MoFEMErrorCode schurSwitchPreconditioner(boost::shared_ptr< BlockStructure > block_mat_data)
Switch preconditioner.
Definition: Schur.cpp:2365
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:215
double
MoFEM::PipelineManager::FaceEle
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Definition: PipelineManager.hpp:35
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
I
constexpr IntegrationType I
Definition: schur_test_diag_mat.cpp:31
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::createOpSchurAssembleBegin
OpSchurAssembleBase * createOpSchurAssembleBegin()
Definition: Schur.cpp:2181
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
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_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
MoFEM::OpBaseImpl::MatSetValuesHook
boost::function< MoFEMErrorCode(ForcesAndSourcesCore::UserDataOperator *op_ptr, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, MatrixDouble &m)> MatSetValuesHook
Definition: FormsIntegrators.hpp:218
MoFEM::DMMoFEMSetBlocMatData
MoFEMErrorCode DMMoFEMSetBlocMatData(DM dm, boost::shared_ptr< BlockStructure >)
Set data for block mat.
Definition: DMMoFEM.cpp:1533
BiLinearForm
SPACE_DIM
constexpr int SPACE_DIM
[Define dimension]
Definition: schur_test_diag_mat.cpp:26
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition: FormsIntegrators.hpp:136
MoFEM::PipelineManager::setDomainLhsIntegrationRule
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Definition: PipelineManager.hpp:503
MoFEM::setSchurA00MatSolvePC
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
Definition: Schur.cpp:2223
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
ElementsAndOps
Definition: child_and_parent.cpp:18
block_mat
SmartPetscObj< Mat > block_mat
Definition: schur_test_diag_mat.cpp:39
DomainEleOp
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition: FormsIntegrators.hpp:104
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
EntData
EntitiesFieldData::EntData EntData
Definition: schur_test_diag_mat.cpp:34
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
debug
constexpr bool debug
Definition: schur_test_diag_mat.cpp:83
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
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
FIELD_DIM
constexpr int FIELD_DIM
Definition: schur_test_diag_mat.cpp:27
main
int main(int argc, char *argv[])
Definition: schur_test_diag_mat.cpp:85
A
constexpr AssemblyType A
Definition: schur_test_diag_mat.cpp:29
MoFEM::BLOCK_SCHUR
@ BLOCK_SCHUR
Definition: FormsIntegrators.hpp:108
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
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::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
MoFEM::SmartPetscObj< Mat >
MoFEM::schurSaveBlockMesh
MoFEMErrorCode schurSaveBlockMesh(boost::shared_ptr< BlockStructure > block_mat_data, std::string filename)
Save block matrix as a mesh.
Definition: Schur.cpp:2380
petsc_mat
SmartPetscObj< Mat > petsc_mat
Definition: schur_test_diag_mat.cpp:38
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
Generated by Doxygen 1.8.17 and hosted at University of Glasgow