v0.14.0
Loading...
Searching...
No Matches
Public Member Functions | Private Types | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
NonlinearPoisson Struct Reference
Collaboration diagram for NonlinearPoisson:
[legend]

Public Member Functions

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

Private Types

using PostProcEle = PostProcBrokenMeshInMoab< FaceEle >
 

Private Member Functions

MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode setIntegrationRules ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
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
 
MPI_Comm mpiComm
 
const int mpiRank
 
SmartPetscObj< DM > dM
 
SmartPetscObj< SNES > snesSolver
 
std::string domainField
 
int order
 
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
 
boost::shared_ptr< FaceEledomainTangentMatrixPipeline
 
boost::shared_ptr< FaceEledomainResidualVectorPipeline
 
boost::shared_ptr< EdgeEleboundaryTangentMatrixPipeline
 
boost::shared_ptr< EdgeEleboundaryResidualVectorPipeline
 
boost::shared_ptr< DataAtGaussPointspreviousUpdate
 
boost::shared_ptr< VectorDouble > fieldValuePtr
 
boost::shared_ptr< MatrixDouble > fieldGradPtr
 
boost::shared_ptr< PostProcElepostProc
 
Range boundaryEntitiesForFieldsplit
 

Detailed Description

Definition at line 10 of file nonlinear_poisson_2d.cpp.

Member Typedef Documentation

◆ PostProcEle

Definition at line 72 of file nonlinear_poisson_2d.cpp.

Constructor & Destructor Documentation

◆ NonlinearPoisson()

NonlinearPoisson::NonlinearPoisson ( MoFEM::Interface m_field)

Definition at line 81 of file nonlinear_poisson_2d.cpp.

82 : domainField("U"), mField(m_field), mpiComm(mField.get_comm()),
84 domainTangentMatrixPipeline = boost::shared_ptr<FaceEle>(new FaceEle(mField));
86 boost::shared_ptr<FaceEle>(new FaceEle(mField));
88 boost::shared_ptr<EdgeEle>(new EdgeEle(mField));
90 boost::shared_ptr<EdgeEle>(new EdgeEle(mField));
91
93 boost::shared_ptr<DataAtGaussPoints>(new DataAtGaussPoints());
94 fieldValuePtr = boost::shared_ptr<VectorDouble>(previousUpdate,
95 &previousUpdate->fieldValue);
96 fieldGradPtr = boost::shared_ptr<MatrixDouble>(previousUpdate,
97 &previousUpdate->fieldGrad);
98}
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
MoFEM::Interface & mField
boost::shared_ptr< EdgeEle > boundaryTangentMatrixPipeline
boost::shared_ptr< VectorDouble > fieldValuePtr
boost::shared_ptr< MatrixDouble > fieldGradPtr
boost::shared_ptr< EdgeEle > boundaryResidualVectorPipeline
boost::shared_ptr< FaceEle > domainResidualVectorPipeline
boost::shared_ptr< DataAtGaussPoints > previousUpdate
boost::shared_ptr< FaceEle > domainTangentMatrixPipeline

Member Function Documentation

◆ assembleSystem()

MoFEMErrorCode NonlinearPoisson::assembleSystem ( )
private

Definition at line 198 of file nonlinear_poisson_2d.cpp.

198 {
200
201 auto det_ptr = boost::make_shared<VectorDouble>();
202 auto jac_ptr = boost::make_shared<MatrixDouble>();
203 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
204
205 { // Push operators to the Pipeline that is responsible for calculating the
206 // domain tangent matrix (LHS)
207
208 // Add default operators to calculate inverse of Jacobian (needed for
209 // implementation of 2D problem but not 3D ones)
210
211 domainTangentMatrixPipeline->getOpPtrVector().push_back(
212 new OpCalculateHOJac<2>(jac_ptr));
213 domainTangentMatrixPipeline->getOpPtrVector().push_back(
214 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
215 domainTangentMatrixPipeline->getOpPtrVector().push_back(
216 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
217 domainTangentMatrixPipeline->getOpPtrVector().push_back(
219
220 // Add default operator to calculate field values at integration points
221 domainTangentMatrixPipeline->getOpPtrVector().push_back(
223 // Add default operator to calculate field gradient at integration points
224 domainTangentMatrixPipeline->getOpPtrVector().push_back(
226
227 // Push operators for domain tangent matrix (LHS)
228 domainTangentMatrixPipeline->getOpPtrVector().push_back(
231 }
232
233 { // Push operators to the Pipeline that is responsible for calculating the
234 // domain residual vector (RHS)
235
236 // Add default operators to calculate inverse of Jacobian (needed for
237 // implementation of 2D problem but not 3D ones)
238 domainResidualVectorPipeline->getOpPtrVector().push_back(
239 new OpCalculateHOJac<2>(jac_ptr));
240 domainResidualVectorPipeline->getOpPtrVector().push_back(
241 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
242 domainResidualVectorPipeline->getOpPtrVector().push_back(
243 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
244 domainResidualVectorPipeline->getOpPtrVector().push_back(
246
247 // Add default operator to calculate field values at integration points
248 domainResidualVectorPipeline->getOpPtrVector().push_back(
250 // Add default operator to calculate field gradient at integration points
251 domainResidualVectorPipeline->getOpPtrVector().push_back(
253
254 domainResidualVectorPipeline->getOpPtrVector().push_back(
257 }
258
259 { // Push operators to the Pipeline that is responsible for calculating the
260 // boundary tangent matrix (LHS)
261
262 boundaryTangentMatrixPipeline->getOpPtrVector().push_back(
264 }
265
266 { // Push operators to the Pipeline that is responsible for calculating
267 // boundary residual vector (RHS)
268
269 // Add default operator to calculate field values at integration points
270 boundaryResidualVectorPipeline->getOpPtrVector().push_back(
272
273 boundaryResidualVectorPipeline->getOpPtrVector().push_back(
276 }
277
278 // get Discrete Manager (SmartPetscObj)
280
281 { // Set operators for nonlinear equations solver (SNES) from MoFEM Pipelines
282
283 // Set operators for calculation of LHS and RHS of domain elements
284 boost::shared_ptr<FaceEle> null_face;
287 null_face);
290 null_face);
291
292 // Set operators for calculation of LHS and RHS of boundary elements
293 boost::shared_ptr<EdgeEle> null_edge;
296 null_edge);
299 null_edge);
300 }
301
303}
@ 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 ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMoFEM.cpp:704
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMoFEM.cpp:745
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition: Simple.hpp:348
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:341
static double boundaryFunction(const double x, const double y, const double z)
static double sourceTermFunction(const double x, const double y, const double z)
SmartPetscObj< DM > dM
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker

◆ boundaryCondition()

MoFEMErrorCode NonlinearPoisson::boundaryCondition ( )
private

Definition at line 157 of file nonlinear_poisson_2d.cpp.

157 {
159
160 // Get boundary edges marked in block named "BOUNDARY_CONDITION"
161 auto get_ents_on_mesh_skin = [&]() {
162 Range boundary_entities;
164 std::string entity_name = it->getName();
165 if (entity_name.compare(0, 18, "BOUNDARY_CONDITION") == 0) {
166 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
167 boundary_entities, true);
168 }
169 }
170 // Add vertices to boundary entities
171 Range boundary_vertices;
172 CHKERR mField.get_moab().get_connectivity(boundary_entities,
173 boundary_vertices, true);
174 boundary_entities.merge(boundary_vertices);
175
176 // Store entities for fieldsplit (block) solver
177 boundaryEntitiesForFieldsplit = boundary_entities;
178
179 return boundary_entities;
180 };
181
182 auto mark_boundary_dofs = [&](Range &&skin_edges) {
183 auto problem_manager = mField.getInterface<ProblemsManager>();
184 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
185 problem_manager->markDofs(simpleInterface->getProblemName(), ROW,
186 skin_edges, *marker_ptr);
187 return marker_ptr;
188 };
189
190 // Get global local vector of marked DOFs. Is global, since is set for all
191 // DOFs on processor. Is local since only DOFs on processor are in the
192 // vector. To access DOFs use local indices.
193 boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
194
196}
@ ROW
Definition: definitions.h:123
@ BLOCKSET
Definition: definitions.h:148
#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:362
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ boundaryFunction()

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

Definition at line 35 of file nonlinear_poisson_2d.cpp.

36 {
37 return sin(x * 10.) * cos(y * 10.);
38 // return 0;
39 }

◆ outputResults()

MoFEMErrorCode NonlinearPoisson::outputResults ( )
private

Definition at line 357 of file nonlinear_poisson_2d.cpp.

357 {
359
360 postProc = boost::make_shared<PostProcEle>(mField);
361
362 auto u_ptr = boost::make_shared<VectorDouble>();
363 postProc->getOpPtrVector().push_back(
365
367
368 postProc->getOpPtrVector().push_back(
369
370 new OpPPMap(postProc->getPostProcMesh(), postProc->getMapGaussPts(),
371
372 {{"U", u_ptr}},
373
374 {}, {}, {}
375
376 )
377
378 );
379
382 boost::dynamic_pointer_cast<FEMethod>(postProc));
383
384 CHKERR postProc->writeFile("out_result.h5m");
385
387}
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:572
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Post post-proc data at points from hash maps.
boost::shared_ptr< PostProcEle > postProc

◆ readMesh()

MoFEMErrorCode NonlinearPoisson::readMesh ( )
private

Definition at line 114 of file nonlinear_poisson_2d.cpp.

114 {
116
120
122}
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180

◆ runProgram()

MoFEMErrorCode NonlinearPoisson::runProgram ( )

Definition at line 100 of file nonlinear_poisson_2d.cpp.

100 {
102
103 readMesh();
104 setupProblem();
108 solveSystem();
110
112}
MoFEMErrorCode outputResults()
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode solveSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode setupProblem()

◆ setIntegrationRules()

MoFEMErrorCode NonlinearPoisson::setIntegrationRules ( )
private

Definition at line 141 of file nonlinear_poisson_2d.cpp.

141 {
143
144 auto domain_rule_lhs = [](int, int, int p) -> int { return 2 * (p + 1); };
145 auto domain_rule_rhs = [](int, int, int p) -> int { return 2 * (p + 1); };
146 domainTangentMatrixPipeline->getRuleHook = domain_rule_lhs;
147 domainResidualVectorPipeline->getRuleHook = domain_rule_rhs;
148
149 auto boundary_rule_lhs = [](int, int, int p) -> int { return 2 * p + 1; };
150 auto boundary_rule_rhs = [](int, int, int p) -> int { return 2 * p + 1; };
151 boundaryTangentMatrixPipeline->getRuleHook = boundary_rule_lhs;
152 boundaryResidualVectorPipeline->getRuleHook = boundary_rule_rhs;
153
155}
static Index< 'p', 3 > p

◆ setupProblem()

MoFEMErrorCode NonlinearPoisson::setupProblem ( )
private

Definition at line 124 of file nonlinear_poisson_2d.cpp.

124 {
126
131
132 int order = 3;
133 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
135
137
139}
@ 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:264
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
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
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611

◆ solveSystem()

MoFEMErrorCode NonlinearPoisson::solveSystem ( )
private

Definition at line 305 of file nonlinear_poisson_2d.cpp.

305 {
307
308 // Create RHS and solution vectors
309 SmartPetscObj<Vec> global_rhs, global_solution;
311 global_solution = vectorDuplicate(global_rhs);
312
313 // Create nonlinear solver (SNES)
315 CHKERR SNESSetFromOptions(snesSolver);
316
317 // Fieldsplit block solver: yes/no
318 if (1) {
319 KSP ksp_solver;
320 CHKERR SNESGetKSP(snesSolver, &ksp_solver);
321 PC pc;
322 CHKERR KSPGetPC(ksp_solver, &pc);
323
324 PetscBool is_pcfs = PETSC_FALSE;
325 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
326
327 // Set up FIELDSPLIT, only when option used -pc_type fieldsplit
328 if (is_pcfs == PETSC_TRUE) {
329 IS is_boundary;
330 const MoFEM::Problem *problem_ptr;
331 CHKERR DMMoFEMGetProblemPtr(dM, &problem_ptr);
332 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
333 problem_ptr->getName(), ROW, domainField, 0, 1, &is_boundary,
335 // CHKERR ISView(is_boundary, PETSC_VIEWER_STDOUT_SELF);
336
337 CHKERR PCFieldSplitSetIS(pc, NULL, is_boundary);
338
339 CHKERR ISDestroy(&is_boundary);
340 }
341 }
342
343 CHKERR SNESSetDM(snesSolver, dM);
344 CHKERR SNESSetUp(snesSolver);
345
346 // Solve the system
347 CHKERR SNESSolve(snesSolver, global_rhs, global_solution);
348 // VecView(global_rhs, PETSC_VIEWER_STDOUT_SELF);
349
350 // Scatter result data on the mesh
351 CHKERR DMoFEMMeshToGlobalVector(dM, global_solution, INSERT_VALUES,
352 SCATTER_REVERSE);
353
355}
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:412
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1153
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:521
auto createSNES(MPI_Comm comm)
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
intrusive_ptr for managing petsc objects
SmartPetscObj< SNES > snesSolver

◆ sourceTermFunction()

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

Definition at line 28 of file nonlinear_poisson_2d.cpp.

29 {
30 return 200 * sin(x * 10.) * cos(y * 10.);
31 // return 1;
32 }

Member Data Documentation

◆ boundaryEntitiesForFieldsplit

Range NonlinearPoisson::boundaryEntitiesForFieldsplit
private

Definition at line 78 of file nonlinear_poisson_2d.cpp.

◆ boundaryMarker

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

Definition at line 59 of file nonlinear_poisson_2d.cpp.

◆ boundaryResidualVectorPipeline

boost::shared_ptr<EdgeEle> NonlinearPoisson::boundaryResidualVectorPipeline
private

Definition at line 65 of file nonlinear_poisson_2d.cpp.

◆ boundaryTangentMatrixPipeline

boost::shared_ptr<EdgeEle> NonlinearPoisson::boundaryTangentMatrixPipeline
private

Definition at line 64 of file nonlinear_poisson_2d.cpp.

◆ dM

SmartPetscObj<DM> NonlinearPoisson::dM
private

Definition at line 51 of file nonlinear_poisson_2d.cpp.

◆ domainField

std::string NonlinearPoisson::domainField
private

Definition at line 55 of file nonlinear_poisson_2d.cpp.

◆ domainResidualVectorPipeline

boost::shared_ptr<FaceEle> NonlinearPoisson::domainResidualVectorPipeline
private

Definition at line 63 of file nonlinear_poisson_2d.cpp.

◆ domainTangentMatrixPipeline

boost::shared_ptr<FaceEle> NonlinearPoisson::domainTangentMatrixPipeline
private

Definition at line 62 of file nonlinear_poisson_2d.cpp.

◆ fieldGradPtr

boost::shared_ptr<MatrixDouble> NonlinearPoisson::fieldGradPtr
private

Definition at line 70 of file nonlinear_poisson_2d.cpp.

◆ fieldValuePtr

boost::shared_ptr<VectorDouble> NonlinearPoisson::fieldValuePtr
private

Definition at line 69 of file nonlinear_poisson_2d.cpp.

◆ mField

MoFEM::Interface& NonlinearPoisson::mField
private

Definition at line 42 of file nonlinear_poisson_2d.cpp.

◆ mpiComm

MPI_Comm NonlinearPoisson::mpiComm
private

Definition at line 46 of file nonlinear_poisson_2d.cpp.

◆ mpiRank

const int NonlinearPoisson::mpiRank
private

Definition at line 48 of file nonlinear_poisson_2d.cpp.

◆ order

int NonlinearPoisson::order
private

Definition at line 56 of file nonlinear_poisson_2d.cpp.

◆ postProc

boost::shared_ptr<PostProcEle> NonlinearPoisson::postProc
private

Definition at line 75 of file nonlinear_poisson_2d.cpp.

◆ previousUpdate

boost::shared_ptr<DataAtGaussPoints> NonlinearPoisson::previousUpdate
private

Definition at line 68 of file nonlinear_poisson_2d.cpp.

◆ simpleInterface

Simple* NonlinearPoisson::simpleInterface
private

Definition at line 43 of file nonlinear_poisson_2d.cpp.

◆ snesSolver

SmartPetscObj<SNES> NonlinearPoisson::snesSolver
private

Definition at line 52 of file nonlinear_poisson_2d.cpp.


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