v0.9.2
lesson3_poisson.cpp
Go to the documentation of this file.
1 /**
2  * \file lesson3_poisson.cpp
3  * \example lesson3_poisson.cpp
4  *
5  * Using PipelineManager interface calculate the divergence of base functions, and
6  * integral of flux on the boundary. Since the h-div space is used, volume
7  * integral and boundary integral should give the same result.
8  */
9 
10 /* This file is part of MoFEM.
11  * MoFEM is free software: you can redistribute it and/or modify it under
12  * the terms of the GNU Lesser General Public License as published by the
13  * Free Software Foundation, either version 3 of the License, or (at your
14  * option) any later version.
15  *
16  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19  * License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23 
24 #include <MoFEM.hpp>
25 
26 using namespace MoFEM;
27 
28 static char help[] = "...\n\n";
29 
30 #include <BasicFiniteElements.hpp>
31 
37 
38 #include <BaseOps.hpp>
39 
44 
45 struct Example {
46 
47  Example(MoFEM::Interface &m_field) : mField(m_field) {}
48 
49  MoFEMErrorCode runProblem();
50 
51 private:
52  MoFEM::Interface &mField;
53 
54  static double approxFunction(const double x, const double y, const double z) {
55  return sin(x * 10.) * cos(y * 10.);
56  }
57  static double nablaFunction(const double x, const double y, const double z) {
58  return 200 * sin(x * 10.) * cos(y * 10.);
59  }
60 
61  static int integrationRule(int, int, int p_data) { return 2 * p_data; };
62 
63  MoFEMErrorCode setUP();
64  MoFEMErrorCode createCommonData();
65  MoFEMErrorCode bC();
66  MoFEMErrorCode OPs();
67  MoFEMErrorCode kspSolve();
68  MoFEMErrorCode postProcess();
69  MoFEMErrorCode checkResults();
70 
72  struct CommonData {
73  boost::shared_ptr<VectorDouble> approxVals;
74  SmartPetscObj<Vec> L2Vec;
75  SmartPetscObj<Vec> resVec;
76  };
77  boost::shared_ptr<CommonData> commonDataPtr;
78  boost::shared_ptr<std::vector<bool>> boundaryMarker;
79 
80  struct OpError : public DomainEleOp {
81  boost::shared_ptr<CommonData> commonDataPtr;
82  OpError(boost::shared_ptr<CommonData> &common_data_ptr)
83  : DomainEleOp("U", OPROW), commonDataPtr(common_data_ptr) {}
84  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
85  };
86 };
87 
90  CHKERR setUP();
91  CHKERR createCommonData();
92  CHKERR bC();
93  CHKERR OPs();
94  CHKERR kspSolve();
95  CHKERR postProcess();
96  CHKERR checkResults();
98 }
99 
100 //! [Set up problem]
103  Simple *simple = mField.getInterface<Simple>();
104  // Add field
105  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
106  CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
107  constexpr int order = 5;
108  CHKERR simple->setFieldOrder("U", order);
109  CHKERR simple->setUp();
111 }
112 //! [Set up problem]
113 
114 //! [Create common data]
117  Simple *simple = mField.getInterface<Simple>();
118  commonDataPtr = boost::make_shared<CommonData>();
119  commonDataPtr->resVec = smartCreateDMVector(simple->getDM());
120  commonDataPtr->L2Vec = createSmartVectorMPI(
121  mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
122  commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
124 }
125 //! [Create common data]
126 
127 //! [Boundary condition]
130 
131  Simple *simple = mField.getInterface<Simple>();
132 
133  auto get_ents_on_mesh_skin = [&]() {
134  Range faces;
135  CHKERR mField.get_moab().get_entities_by_type(0, MBTRI, faces);
136  Skinner skin(&mField.get_moab());
137  Range skin_edges;
138  CHKERR skin.find_skin(0, faces, false, skin_edges);
139  Range skin_verts;
140  CHKERR mField.get_moab().get_connectivity(skin_edges, skin_verts, true);
141  skin_edges.merge(skin_verts);
142  return skin_edges;
143  };
144 
145  auto mark_boundary_dofs = [&](Range &&skin_edges) {
146  auto problem_manager = mField.getInterface<ProblemsManager>();
147  auto marker_ptr = boost::make_shared<std::vector<bool>>();
148  problem_manager->markDofs(simple->getProblemName(), ROW, skin_edges,
149  *marker_ptr);
150  return marker_ptr;
151  };
152 
153  // Get global local vector of marked DOFs. Is global, since is set for all
154  // DOFs on processor. Is local since only DOFs on processor are in the
155  // vector. To access DOFs use local indices.
156  boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
157 
159 }
160 //! [Boundary condition]
161 
162 //! [Push operators to pipeline]
165  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
166  pipeline_mng->getOpDomainLhsPipeline().push_back(
168  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetInvJacH1ForFace(invJac));
169  auto beta = [](const double, const double, const double) { return 1; };
170  pipeline_mng->getOpDomainLhsPipeline().push_back(
171  new OpDomainGradGrad("U", "U", beta, boundaryMarker));
172  pipeline_mng->getOpDomainRhsPipeline().push_back(
173  new OpDomainSource("U", Example::nablaFunction, boundaryMarker));
174  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integrationRule);
175  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integrationRule);
176  pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpBoundaryMass("U", "U", beta));
177  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
178  new OpBoundarySource("U", approxFunction));
179  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integrationRule);
180  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integrationRule);
182 }
183 //! [Push operators to pipeline]
184 
185 //! [Solve]
188  Simple *simple = mField.getInterface<Simple>();
189  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
190  auto solver = pipeline_mng->createKSP();
191  CHKERR KSPSetFromOptions(solver);
192  CHKERR KSPSetUp(solver);
193 
194  auto dm = simple->getDM();
195  auto D = smartCreateDMVector(dm);
196  auto F = smartVectorDuplicate(D);
197 
198  CHKERR KSPSolve(solver, F, D);
199  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
200  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
201  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
203 }
204 
205 //! [Solve]
208  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
209  pipeline_mng->getDomainLhsFE().reset();
210  pipeline_mng->getBoundaryLhsFE().reset();
211  pipeline_mng->getBoundaryRhsFE().reset();
212  auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
213  post_proc_fe->generateReferenceElementMesh();
214  post_proc_fe->addFieldValuesPostProc("U");
215  pipeline_mng->getDomainRhsFE() = post_proc_fe;
216  CHKERR pipeline_mng->loopFiniteElements();
217  CHKERR post_proc_fe->writeFile("out_poisson.h5m");
219 }
220 //! [Postprocess results]
221 
222 //! [Solve]
225  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
226  pipeline_mng->getDomainLhsFE().reset();
227  pipeline_mng->getDomainRhsFE().reset();
228  pipeline_mng->getBoundaryLhsFE().reset();
229  pipeline_mng->getBoundaryRhsFE().reset();
230  pipeline_mng->getOpDomainRhsPipeline().clear();
231  pipeline_mng->getOpDomainRhsPipeline().push_back(
232  new OpCalculateScalarFieldValues("U", commonDataPtr->approxVals));
233  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpError(commonDataPtr));
234  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integrationRule);
235  CHKERR pipeline_mng->loopFiniteElements();
236  CHKERR VecAssemblyBegin(commonDataPtr->L2Vec);
237  CHKERR VecAssemblyEnd(commonDataPtr->L2Vec);
238  CHKERR VecAssemblyBegin(commonDataPtr->resVec);
239  CHKERR VecAssemblyEnd(commonDataPtr->resVec);
240  double nrm2;
241  CHKERR VecNorm(commonDataPtr->resVec, NORM_2, &nrm2);
242  const double *array;
243  CHKERR VecGetArrayRead(commonDataPtr->L2Vec, &array);
244  if (mField.get_comm_rank() == 0)
245  PetscPrintf(PETSC_COMM_SELF, "Error %6.4e Vec norm %6.4e\n", sqrt(array[0]),
246  nrm2);
247  CHKERR VecRestoreArrayRead(commonDataPtr->L2Vec, &array);
248  constexpr double eps = 1e-8;
249  if (nrm2 > eps)
250  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
251  "Not converged solution");
253 }
254 //! [Solver]
255 
256 int main(int argc, char *argv[]) {
257 
258  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
259 
260  try {
261 
262  //! [Register MoFEM discrete manager in PETSc]
263  DMType dm_name = "DMMOFEM";
264  CHKERR DMRegister_MoFEM(dm_name);
265  //! [Register MoFEM discrete manager in PETSc
266 
267  //! [Create MoAB]
268  moab::Core mb_instance; ///< mesh database
269  moab::Interface &moab = mb_instance; ///< mesh database interface
270  //! [Create MoAB]
271 
272  //! [Create MoFEM]
273  MoFEM::Core core(moab); ///< finite element database
274  MoFEM::Interface &m_field = core; ///< finite element database insterface
275  //! [Create MoFEM]
276 
277  //! [Load mesh]
278  Simple *simple = m_field.getInterface<Simple>();
279  CHKERR simple->getOptions();
280  CHKERR simple->loadFile("");
281  //! [Load mesh]
282 
283  //! [Example]
284  Example ex(m_field);
285  CHKERR ex.runProblem();
286  //! [Example]
287  }
288  CATCH_ERRORS;
289 
291 }
292 
293 MoFEMErrorCode Example::OpError::doWork(int side, EntityType type,
294  EntData &data) {
296 
297  if (const size_t nb_dofs = data.getIndices().size()) {
298 
299  const int nb_integration_pts = getGaussPts().size2();
300  auto t_w = getFTensor0IntegrationWeight();
301  auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
302  auto t_coords = getFTensor1CoordsAtGaussPts();
303 
304  VectorDouble nf(nb_dofs, false);
305  nf.clear();
306 
308  const double volume = getMeasure();
309 
310  auto t_row_base = data.getFTensor0N();
311  double error = 0;
312  for (int gg = 0; gg != nb_integration_pts; ++gg) {
313 
314  const double alpha = t_w * volume;
315  double diff = t_val - Example::approxFunction(t_coords(0), t_coords(1),
316  t_coords(2));
317  error += alpha * pow(diff, 2);
318 
319  for (size_t r = 0; r != nb_dofs; ++r) {
320  nf[r] += alpha * t_row_base * diff;
321  ++t_row_base;
322  }
323 
324  ++t_w;
325  ++t_val;
326  ++t_coords;
327  }
328 
329  const int index = 0;
330  CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
331  CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
332  }
333 
335 }
OpTools< DomainEleOp >::OpSource< 2 > OpDomainSource
OpTools< BoundaryEleOp >::OpSource< 2 > OpBoundarySource
boost::shared_ptr< FEMethod > & getDomainLhsFE()
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:141
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Deprecated interface functions.
Problem manager is used to build and partition problems.
MoFEMErrorCode createCommonData()
[Set up problem]
MoFEMErrorCode kspSolve()
[Push operators to pipeline]
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
static double approxFunction(const double x, const double y, const double z)
[Grad residual operator]
Definition: BaseOps.hpp:239
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:445
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MatrixDouble invJac
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
Example(MoFEM::Interface &m_field)
Core (interface) class.
Definition: Core.hpp:50
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
Definition: AuxPETSc.hpp:238
OpError(boost::shared_ptr< CommonData > &common_data_ptr)
OpTools< BoundaryEleOp >::OpMass OpBoundaryMass
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
static char help[]
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::shared_ptr< std::vector< bool > > boundaryMarker
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:939
auto createSmartVectorMPI
Create MPI Vector.
Definition: AuxPETSc.hpp:224
const VectorInt & getIndices() const
Get global indices of dofs on entity.
MoFEMErrorCode runProblem()
[Operator]
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEMErrorCode OPs()
[Boundary condition]
int main(int argc, char *argv[])
[Solver]
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
MoFEMErrorCode checkResults()
[Print results]
Calculate inverse of jacobian for face element.
MoFEMErrorCode bC()
[Create common data]
[Source operator]
Definition: BaseOps.hpp:97
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
static moab::Error error
static int integrationRule(int, int, int p_data)
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
PipelineManager interface.
#define CHKERR
Inline error check.
Definition: definitions.h:602
boost::shared_ptr< FEMethod > & getBoundaryLhsFE()
MoFEMErrorCode setUP()
[Set up problem]
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
static double nablaFunction(const double x, const double y, const double z)
Transform local reference derivatives of shape functions to global derivatives.
OpTools< DomainEleOp >::OpGradGrad< 2 > OpDomainGradGrad
constexpr int order
MoFEMErrorCode postProcess()
[Integrate]
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
boost::ptr_vector< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
DataForcesAndSourcesCore::EntData EntData
[Source operator]
Definition: BaseOps.hpp:137
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const Range ents, std::vector< bool > &marker)
Create vector with marked indices.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:439
continuous field
Definition: definitions.h:177
static const double eps
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
Get value at integration points for scalar field.
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
MatrixDouble invJac
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:69
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26