v0.14.0
poisson_2d_lagrange_multiplier.cpp
Go to the documentation of this file.
1 #include <stdlib.h>
4 
5 using namespace MoFEM;
7 
8 using PostProcFaceEle =
10 using PostProcEdgeEle =
12 
13 static char help[] = "...\n\n";
14 
16 public:
18 
19  // Declaration of the main function to run analysis
20  MoFEMErrorCode runProgram();
21 
22 private:
23  // Declaration of other main functions called in runProgram()
24  MoFEMErrorCode readMesh();
25  MoFEMErrorCode setupProblem();
26  MoFEMErrorCode boundaryCondition();
27  MoFEMErrorCode assembleSystem();
28  MoFEMErrorCode setIntegrationRules();
29  MoFEMErrorCode solveSystem();
30  MoFEMErrorCode outputResults();
31 
32  // Function to calculate the Source term
33  static double sourceTermFunction(const double x, const double y,
34  const double z) {
35  return 200 * sin(x * 10.) * cos(y * 10.);
36  // return 1;
37  }
38  // Function to calculate the Boundary term
39  static double boundaryFunction(const double x, const double y,
40  const double z) {
41  return sin(x * 10.) * cos(y * 10.);
42  // return 0;
43  }
44 
45  // Main interfaces
48 
49  // Field name and approximation order
50  std::string domainField; // displacement field
51  std::string boundaryField; // Lagrange multiplier field
52  int oRder;
53 
54  // Object to mark boundary entities for the assembling of domain elements
55  boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
56 
57  // Boundary entities marked for fieldsplit (block) solver - optional
59 };
60 
62  MoFEM::Interface &m_field)
63  : domainField("U"), boundaryField("L"), mField(m_field) {}
64 
67 
71 
73 }
74 
77 
82 
83  int oRder = 3;
84  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
87 
89 
91 }
92 
95 
96  // Get boundary edges marked in block named "BOUNDARY_CONDITION"
97  auto get_ents_on_mesh = [&]() {
98  Range boundary_entities;
100  std::string entity_name = it->getName();
101  if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
102  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
103  boundary_entities, true);
104  }
105  }
106  // Add vertices to boundary entities
107  Range boundary_vertices;
108  CHKERR mField.get_moab().get_connectivity(boundary_entities,
109  boundary_vertices, true);
110  boundary_entities.merge(boundary_vertices);
111 
112  // Store entities for fieldsplit (block) solver
113  boundaryEntitiesForFieldsplit = boundary_entities;
114 
115  return boundary_entities;
116  };
117 
118  auto mark_boundary_dofs = [&](Range &&skin_edges) {
119  auto problem_manager = mField.getInterface<ProblemsManager>();
120  auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
121  problem_manager->markDofs(simpleInterface->getProblemName(), ROW,
122  skin_edges, *marker_ptr);
123  return marker_ptr;
124  };
125 
126  // Get global local vector of marked DOFs. Is global, since is set for all
127  // DOFs on processor. Is local since only DOFs on processor are in the
128  // vector. To access DOFs use local indices.
129  boundaryMarker = mark_boundary_dofs(get_ents_on_mesh());
130 
132 }
133 
136 
137  auto pipeline_mng = mField.getInterface<PipelineManager>();
138  auto det_ptr = boost::make_shared<VectorDouble>();
139  auto jac_ptr = boost::make_shared<MatrixDouble>();
140  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
141 
142  { // Push operators to the Pipeline that is responsible for calculating LHS of
143  // domain elements
144  pipeline_mng->getOpDomainLhsPipeline().push_back(
145  new OpCalculateHOJac<2>(jac_ptr));
146  pipeline_mng->getOpDomainLhsPipeline().push_back(
147  new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
148  pipeline_mng->getOpDomainLhsPipeline().push_back(
149  new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
150  pipeline_mng->getOpDomainLhsPipeline().push_back(
151  new OpSetHOWeightsOnFace());
152  pipeline_mng->getOpDomainLhsPipeline().push_back(
154  }
155 
156  { // Push operators to the Pipeline that is responsible for calculating RHS of
157  // domain elements
158  pipeline_mng->getOpDomainRhsPipeline().push_back(
160  }
161 
162  { // Push operators to the Pipeline that is responsible for calculating LHS of
163  // boundary elements (Lagrange multiplier)
164  pipeline_mng->getOpBoundaryLhsPipeline().push_back(
166  pipeline_mng->getOpBoundaryLhsPipeline().push_back(
168  }
169 
170  { // Push operators to the Pipeline that is responsible for calculating RHS of
171  // boundary elements (Lagrange multiplier)
172  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
174  }
175 
177 }
178 
181 
182  auto pipeline_mng = mField.getInterface<PipelineManager>();
183 
184  auto domain_rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
185  auto domain_rule_rhs = [](int, int, int p) -> int { return 2 * (p - 1); };
186  CHKERR pipeline_mng->setDomainLhsIntegrationRule(domain_rule_lhs);
187  CHKERR pipeline_mng->setDomainRhsIntegrationRule(domain_rule_rhs);
188 
189  auto boundary_rule_lhs = [](int, int, int p) -> int { return 2 * p; };
190  auto boundary_rule_rhs = [](int, int, int p) -> int { return 2 * p; };
191  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_lhs);
192  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_rhs);
193 
195 }
196 
199 
200  auto pipeline_mng = mField.getInterface<PipelineManager>();
201 
202  auto ksp_solver = pipeline_mng->createKSP();
203  CHKERR KSPSetFromOptions(ksp_solver);
204 
205  // Create RHS and solution vectors
206  auto dm = simpleInterface->getDM();
207  auto F = createDMVector(dm);
208  auto D = vectorDuplicate(F);
209 
210  // Setup fieldsplit (block) solver - optional: yes/no
211  if (1) {
212  PC pc;
213  CHKERR KSPGetPC(ksp_solver, &pc);
214  PetscBool is_pcfs = PETSC_FALSE;
215  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
216 
217  // Set up FIELDSPLIT, only when user set -pc_type fieldsplit
218  // Identify the index for boundary entities, remaining will be for domain
219  // Then split the fields for boundary and domain for solving
220  if (is_pcfs == PETSC_TRUE) {
221  IS is_domain, is_boundary;
222  cerr << "Running FIELDSPLIT..." << endl;
223  const MoFEM::Problem *problem_ptr;
224  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
225  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
226  problem_ptr->getName(), ROW, domainField, 0, 1, &is_boundary,
228  // CHKERR ISView(is_boundary, PETSC_VIEWER_STDOUT_SELF);
229 
230  CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
231 
232  CHKERR ISDestroy(&is_boundary);
233  }
234  }
235 
236  CHKERR KSPSetUp(ksp_solver);
237 
238  // Solve the system
239  CHKERR KSPSolve(ksp_solver, F, D);
240 
241  // Scatter result data on the mesh
242  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
243  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
244  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
245 
247 }
248 
251 
252  auto pipeline_mng = mField.getInterface<PipelineManager>();
253  pipeline_mng->getDomainLhsFE().reset();
254  pipeline_mng->getBoundaryLhsFE().reset();
255 
256  auto d_ptr = boost::make_shared<VectorDouble>();
257  auto l_ptr = boost::make_shared<VectorDouble>();
258 
260 
261  auto post_proc_domain_fe = boost::make_shared<PostProcFaceEle>(mField);
262  post_proc_domain_fe->getOpPtrVector().push_back(
264  post_proc_domain_fe->getOpPtrVector().push_back(
265  new OpPPMap(post_proc_domain_fe->getPostProcMesh(),
266  post_proc_domain_fe->getMapGaussPts(), {{domainField, d_ptr}},
267  {}, {}, {}));
268  pipeline_mng->getDomainRhsFE() = post_proc_domain_fe;
269 
270  auto post_proc_boundary_fe = boost::make_shared<PostProcEdgeEle>(mField);
271  post_proc_boundary_fe->getOpPtrVector().push_back(
272  new OpCalculateScalarFieldValues(boundaryField, l_ptr));
273  post_proc_boundary_fe->getOpPtrVector().push_back(
274  new OpPPMap(post_proc_boundary_fe->getPostProcMesh(),
275  post_proc_boundary_fe->getMapGaussPts(),
276  {{boundaryField, l_ptr}}, {}, {}, {}));
277  pipeline_mng->getBoundaryRhsFE() = post_proc_boundary_fe;
278 
279  CHKERR pipeline_mng->loopFiniteElements();
280  CHKERR post_proc_domain_fe->writeFile("out_result_domain.h5m");
281  CHKERR post_proc_boundary_fe->writeFile("out_result_boundary.h5m");
282 
284 }
285 
288 
289  readMesh();
290  setupProblem();
292  assembleSystem();
294  solveSystem();
295  outputResults();
296 
298 }
299 
300 int main(int argc, char *argv[]) {
301 
302  // Initialisation of MoFEM/PETSc and MOAB data structures
303  const char param_file[] = "param_file.petsc";
304  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
305 
306  // Error handling
307  try {
308  // Register MoFEM discrete manager in PETSc
309  DMType dm_name = "DMMOFEM";
310  CHKERR DMRegister_MoFEM(dm_name);
311 
312  // Create MOAB instance
313  moab::Core mb_instance; // mesh database
314  moab::Interface &moab = mb_instance; // mesh database interface
315 
316  // Create MoFEM instance
317  MoFEM::Core core(moab); // finite element database
318  MoFEM::Interface &m_field = core; // finite element interface
319 
320  // Run the main analysis
321  Poisson2DLagrangeMultiplier poisson_problem(m_field);
322  CHKERR poisson_problem.runProgram();
323  }
324  CATCH_ERRORS;
325 
326  // Finish work: cleaning memory, getting statistics, etc.
328 
329  return 0;
330 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
Poisson2DLagrangeMultiplier::sourceTermFunction
static double sourceTermFunction(const double x, const double y, const double z)
Definition: poisson_2d_lagrange_multiplier.cpp:33
MoFEM::PostProcBrokenMeshInMoab< EdgeElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:684
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
Poisson2DLagrangeMultiplier::assembleSystem
MoFEMErrorCode assembleSystem()
Definition: poisson_2d_lagrange_multiplier.cpp:134
Poisson2DLagrangeMultiplier::boundaryField
std::string boundaryField
Definition: poisson_2d_lagrange_multiplier.cpp:51
Poisson2DLagrangeMultiplier::solveSystem
MoFEMErrorCode solveSystem()
Definition: poisson_2d_lagrange_multiplier.cpp:197
main
int main(int argc, char *argv[])
Definition: poisson_2d_lagrange_multiplier.cpp:300
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpSetHOInvJacToScalarBases< 2 >
Definition: HODataOperators.hpp:78
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
BasicFiniteElements.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::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::OpCalculateHOJac< 2 >
Definition: HODataOperators.hpp:273
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:123
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
Poisson2DLagrangeMultiplier::Poisson2DLagrangeMultiplier
Poisson2DLagrangeMultiplier(MoFEM::Interface &m_field)
Definition: poisson_2d_lagrange_multiplier.cpp:61
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1067
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
Poisson2DLagrangeMultiplier::readMesh
MoFEMErrorCode readMesh()
Definition: poisson_2d_lagrange_multiplier.cpp:65
Poisson2DLagrangeMultiplier::outputResults
MoFEMErrorCode outputResults()
Definition: poisson_2d_lagrange_multiplier.cpp:249
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
Poisson2DLagrangeMultiplier::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
Definition: poisson_2d_lagrange_multiplier.cpp:179
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Poisson2DLagrangeMultiplierOperators::OpBoundaryLhsC
Definition: poisson_2d_lagrange_multiplier.hpp:169
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
Poisson2DLagrangeMultiplier::boundaryEntitiesForFieldsplit
Range boundaryEntitiesForFieldsplit
Definition: poisson_2d_lagrange_multiplier.cpp:58
Poisson2DLagrangeMultiplier::boundaryFunction
static double boundaryFunction(const double x, const double y, const double z)
Definition: poisson_2d_lagrange_multiplier.cpp:39
Poisson2DLagrangeMultiplier::boundaryCondition
MoFEMErrorCode boundaryCondition()
Definition: poisson_2d_lagrange_multiplier.cpp:93
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::Simple::addDomainField
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:264
MoFEM::PipelineManager::createKSP
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
Definition: PipelineManager.cpp:49
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
Poisson2DLagrangeMultiplier::simpleInterface
Simple * simpleInterface
Definition: poisson_2d_lagrange_multiplier.cpp:47
MoFEM::OpSetHOWeightsOnFace
Modify integration weights on face to take in account higher-order geometry.
Definition: HODataOperators.hpp:122
MoFEM::ProblemsManager::markDofs
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const enum MarkOP op, const Range ents, std::vector< unsigned char > &marker) const
Create vector with marked indices.
Definition: ProblemsManager.cpp:3200
help
static char help[]
Definition: poisson_2d_lagrange_multiplier.cpp:13
Poisson2DLagrangeMultiplier::setupProblem
MoFEMErrorCode setupProblem()
Definition: poisson_2d_lagrange_multiplier.cpp:75
Poisson2DLagrangeMultiplier::domainField
std::string domainField
Definition: poisson_2d_lagrange_multiplier.cpp:50
Poisson2DLagrangeMultiplierOperators::OpDomainRhsF
Definition: poisson_2d_lagrange_multiplier.hpp:104
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
Poisson2DLagrangeMultiplierOperators
Definition: poisson_2d_lagrange_multiplier.hpp:15
poisson_2d_lagrange_multiplier.hpp
Range
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401
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::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
Poisson2DLagrangeMultiplier
Definition: poisson_2d_lagrange_multiplier.cpp:15
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:426
Poisson2DLagrangeMultiplier::runProgram
MoFEMErrorCode runProgram()
Definition: poisson_2d_lagrange_multiplier.cpp:286
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3254
MoFEM::Simple::addBoundaryField
MoFEMErrorCode addBoundaryField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on boundary.
Definition: Simple.cpp:282
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
Poisson2DLagrangeMultiplier::oRder
int oRder
Definition: poisson_2d_lagrange_multiplier.cpp:52
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
MoFEM::Problem::getName
auto getName() const
Definition: ProblemsMultiIndices.hpp:372
Poisson2DLagrangeMultiplier::mField
MoFEM::Interface & mField
Definition: poisson_2d_lagrange_multiplier.cpp:46
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:362
convert.int
int
Definition: convert.py:64
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:416
Poisson2DLagrangeMultiplierOperators::OpBoundaryRhsG
Definition: poisson_2d_lagrange_multiplier.hpp:236
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
Poisson2DLagrangeMultiplierOperators::OpDomainLhsK
Definition: poisson_2d_lagrange_multiplier.hpp:23
Poisson2DLagrangeMultiplier::boundaryMarker
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: poisson_2d_lagrange_multiplier.cpp:55
F
@ F
Definition: free_surface.cpp:394
MoFEM::PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:677
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698