v0.15.0
Loading...
Searching...
No Matches
Poisson2DLagrangeMultiplier Struct Reference
Collaboration diagram for Poisson2DLagrangeMultiplier:
[legend]

Public Member Functions

 Poisson2DLagrangeMultiplier (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProgram ()
 

Private Member Functions

MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode setIntegrationRules ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 

Static Private Member Functions

static double sourceTermFunction (const double x, const double y, const double z)
 
static double boundaryFunction (const double x, const double y, const double z)
 

Private Attributes

MoFEM::InterfacemField
 
SimplesimpleInterface
 
std::string domainField
 
std::string boundaryField
 
int oRder
 
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
 
Range boundaryEntitiesForFieldsplit
 

Detailed Description

Definition at line 16 of file poisson_2d_lagrange_multiplier.cpp.

Constructor & Destructor Documentation

◆ Poisson2DLagrangeMultiplier()

Poisson2DLagrangeMultiplier::Poisson2DLagrangeMultiplier ( MoFEM::Interface & m_field)

Member Function Documentation

◆ assembleSystem()

MoFEMErrorCode Poisson2DLagrangeMultiplier::assembleSystem ( )
private

Definition at line 135 of file poisson_2d_lagrange_multiplier.cpp.

135 {
137
138 auto pipeline_mng = mField.getInterface<PipelineManager>();
139 auto det_ptr = boost::make_shared<VectorDouble>();
140 auto jac_ptr = boost::make_shared<MatrixDouble>();
141 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
142
143 { // Push operators to the Pipeline that is responsible for calculating LHS of
144 // domain elements
145 pipeline_mng->getOpDomainLhsPipeline().push_back(
146 new OpCalculateHOJac<2>(jac_ptr));
147 pipeline_mng->getOpDomainLhsPipeline().push_back(
148 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
149 pipeline_mng->getOpDomainLhsPipeline().push_back(
150 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
151 pipeline_mng->getOpDomainLhsPipeline().push_back(
153 pipeline_mng->getOpDomainLhsPipeline().push_back(
155 }
156
157 { // Push operators to the Pipeline that is responsible for calculating RHS of
158 // domain elements
159 pipeline_mng->getOpDomainRhsPipeline().push_back(
161 }
162
163 { // Push operators to the Pipeline that is responsible for calculating LHS of
164 // boundary elements (Lagrange multiplier)
165 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
167 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
169 }
170
171 { // Push operators to the Pipeline that is responsible for calculating RHS of
172 // boundary elements (Lagrange multiplier)
173 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
175 }
176
178}
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
static double sourceTermFunction(const double x, const double y, const double z)
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
static double boundaryFunction(const double x, const double y, const double z)

◆ boundaryCondition()

MoFEMErrorCode Poisson2DLagrangeMultiplier::boundaryCondition ( )
private

Definition at line 94 of file poisson_2d_lagrange_multiplier.cpp.

94 {
96
97 // Get boundary edges marked in block named "BOUNDARY_CONDITION"
98 auto get_ents_on_mesh = [&]() {
99 Range boundary_entities;
101 std::string entity_name = it->getName();
102 if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
103 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
104 boundary_entities, true);
105 }
106 }
107 // Add vertices to boundary entities
108 Range boundary_vertices;
109 CHKERR mField.get_moab().get_connectivity(boundary_entities,
110 boundary_vertices, true);
111 boundary_entities.merge(boundary_vertices);
112
113 // Store entities for fieldsplit (block) solver
114 boundaryEntitiesForFieldsplit = boundary_entities;
115
116 return boundary_entities;
117 };
118
119 auto mark_boundary_dofs = [&](Range &&skin_edges) {
120 auto problem_manager = mField.getInterface<ProblemsManager>();
121 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
122 problem_manager->markDofs(simpleInterface->getProblemName(), ROW,
123 skin_edges, *marker_ptr);
124 return marker_ptr;
125 };
126
127 // Get global local vector of marked DOFs. Is global, since is set for all
128 // DOFs on processor. Is local since only DOFs on processor are in the
129 // vector. To access DOFs use local indices.
130 boundaryMarker = mark_boundary_dofs(get_ents_on_mesh());
131
133}
@ ROW
@ BLOCKSET
#define CHKERR
Inline error check.
#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.
virtual moab::Interface & get_moab()=0
Problem manager is used to build and partition problems.
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.
const std::string getProblemName() const
Get the Problem Name.
Definition Simple.hpp:410

◆ boundaryFunction()

static double Poisson2DLagrangeMultiplier::boundaryFunction ( const double x,
const double y,
const double z )
inlinestaticprivate

Definition at line 40 of file poisson_2d_lagrange_multiplier.cpp.

41 {
42 return sin(x * 10.) * cos(y * 10.);
43 // return 0;
44 }

◆ outputResults()

MoFEMErrorCode Poisson2DLagrangeMultiplier::outputResults ( )
private

Definition at line 250 of file poisson_2d_lagrange_multiplier.cpp.

250 {
252
253 auto pipeline_mng = mField.getInterface<PipelineManager>();
254 pipeline_mng->getDomainLhsFE().reset();
255 pipeline_mng->getBoundaryLhsFE().reset();
256
257 auto d_ptr = boost::make_shared<VectorDouble>();
258 auto l_ptr = boost::make_shared<VectorDouble>();
259
261
262 auto post_proc_domain_fe = boost::make_shared<PostProcFaceEle>(mField);
263 post_proc_domain_fe->getOpPtrVector().push_back(
265 post_proc_domain_fe->getOpPtrVector().push_back(
266 new OpPPMap(post_proc_domain_fe->getPostProcMesh(),
267 post_proc_domain_fe->getMapGaussPts(), {{domainField, d_ptr}},
268 {}, {}, {}));
269 pipeline_mng->getDomainRhsFE() = post_proc_domain_fe;
270
271 auto post_proc_boundary_fe = boost::make_shared<PostProcEdgeEle>(mField);
272 post_proc_boundary_fe->getOpPtrVector().push_back(
274 post_proc_boundary_fe->getOpPtrVector().push_back(
275 new OpPPMap(post_proc_boundary_fe->getPostProcMesh(),
276 post_proc_boundary_fe->getMapGaussPts(),
277 {{boundaryField, l_ptr}}, {}, {}, {}));
278 pipeline_mng->getBoundaryRhsFE() = post_proc_boundary_fe;
279
280 CHKERR pipeline_mng->loopFiniteElements();
281 CHKERR post_proc_domain_fe->writeFile("out_result_domain.h5m");
282 CHKERR post_proc_boundary_fe->writeFile("out_result_boundary.h5m");
283
285}
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
boost::shared_ptr< FEMethod > & getDomainLhsFE()

◆ readMesh()

MoFEMErrorCode Poisson2DLagrangeMultiplier::readMesh ( )
private

Definition at line 66 of file poisson_2d_lagrange_multiplier.cpp.

66 {
68
72
74}
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:191
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180

◆ runProgram()

◆ setIntegrationRules()

MoFEMErrorCode Poisson2DLagrangeMultiplier::setIntegrationRules ( )
private

Definition at line 180 of file poisson_2d_lagrange_multiplier.cpp.

180 {
182
183 auto pipeline_mng = mField.getInterface<PipelineManager>();
184
185 auto domain_rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
186 auto domain_rule_rhs = [](int, int, int p) -> int { return 2 * (p - 1); };
187 CHKERR pipeline_mng->setDomainLhsIntegrationRule(domain_rule_lhs);
188 CHKERR pipeline_mng->setDomainRhsIntegrationRule(domain_rule_rhs);
189
190 auto boundary_rule_lhs = [](int, int, int p) -> int { return 2 * p; };
191 auto boundary_rule_rhs = [](int, int, int p) -> int { return 2 * p; };
192 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_lhs);
193 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(boundary_rule_rhs);
194
196}
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)

◆ setupProblem()

MoFEMErrorCode Poisson2DLagrangeMultiplier::setupProblem ( )
private

Definition at line 76 of file poisson_2d_lagrange_multiplier.cpp.

76 {
78
83
84 int oRder = 3;
85 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &oRder, PETSC_NULLPTR);
88
90
92}
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
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:261
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:355
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:575
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:736

◆ solveSystem()

MoFEMErrorCode Poisson2DLagrangeMultiplier::solveSystem ( )
private

Definition at line 198 of file poisson_2d_lagrange_multiplier.cpp.

198 {
200
201 auto pipeline_mng = mField.getInterface<PipelineManager>();
202
203 auto ksp_solver = pipeline_mng->createKSP();
204 CHKERR KSPSetFromOptions(ksp_solver);
205
206 // Create RHS and solution vectors
207 auto dm = simpleInterface->getDM();
208 auto F = createDMVector(dm);
209 auto D = vectorDuplicate(F);
210
211 // Setup fieldsplit (block) solver - optional: yes/no
212 if (1) {
213 PC pc;
214 CHKERR KSPGetPC(ksp_solver, &pc);
215 PetscBool is_pcfs = PETSC_FALSE;
216 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
217
218 // Set up FIELDSPLIT, only when user set -pc_type fieldsplit
219 // Identify the index for boundary entities, remaining will be for domain
220 // Then split the fields for boundary and domain for solving
221 if (is_pcfs == PETSC_TRUE) {
222 IS is_boundary;
223 cerr << "Running FIELDSPLIT..." << endl;
224 const MoFEM::Problem *problem_ptr;
225 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
226 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
227 problem_ptr->getName(), ROW, domainField, 0, 1, &is_boundary,
229 // CHKERR ISView(is_boundary, PETSC_VIEWER_STDOUT_SELF);
230
231 CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
232
233 CHKERR ISDestroy(&is_boundary);
234 }
235 }
236
237 CHKERR KSPSetUp(ksp_solver);
238
239 // Solve the system
240 CHKERR KSPSolve(ksp_solver, F, D);
241
242 // Scatter result data on the mesh
243 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
244 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
245 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
246
248}
@ F
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:514
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition DMMoFEM.cpp:422
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
double D
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
keeps basic data about problem
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:800

◆ sourceTermFunction()

static double Poisson2DLagrangeMultiplier::sourceTermFunction ( const double x,
const double y,
const double z )
inlinestaticprivate

Definition at line 34 of file poisson_2d_lagrange_multiplier.cpp.

35 {
36 return 200 * sin(x * 10.) * cos(y * 10.);
37 // return 1;
38 }

Member Data Documentation

◆ boundaryEntitiesForFieldsplit

Range Poisson2DLagrangeMultiplier::boundaryEntitiesForFieldsplit
private

Definition at line 59 of file poisson_2d_lagrange_multiplier.cpp.

◆ boundaryField

std::string Poisson2DLagrangeMultiplier::boundaryField
private

Definition at line 52 of file poisson_2d_lagrange_multiplier.cpp.

◆ boundaryMarker

boost::shared_ptr<std::vector<unsigned char> > Poisson2DLagrangeMultiplier::boundaryMarker
private

Definition at line 56 of file poisson_2d_lagrange_multiplier.cpp.

◆ domainField

std::string Poisson2DLagrangeMultiplier::domainField
private

Definition at line 51 of file poisson_2d_lagrange_multiplier.cpp.

◆ mField

MoFEM::Interface& Poisson2DLagrangeMultiplier::mField
private

Definition at line 47 of file poisson_2d_lagrange_multiplier.cpp.

◆ oRder

int Poisson2DLagrangeMultiplier::oRder
private

Definition at line 53 of file poisson_2d_lagrange_multiplier.cpp.

◆ simpleInterface

Simple* Poisson2DLagrangeMultiplier::simpleInterface
private

Definition at line 48 of file poisson_2d_lagrange_multiplier.cpp.


The documentation for this struct was generated from the following file: