v0.10.0
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 // TODO: [CORE-62] Improve method of enforcing boundary condition
25 
26 #include <MoFEM.hpp>
27 
28 using namespace MoFEM;
29 
30 static char help[] = "...\n\n";
31 
32 #include <BasicFiniteElements.hpp>
33 
39 
40 #include <BaseOps.hpp>
41 
46 
47 struct Example {
48 
49  Example(MoFEM::Interface &m_field) : mField(m_field) {}
50 
51  MoFEMErrorCode runProblem();
52 
53 private:
55 
56  static double approxFunction(const double x, const double y, const double z) {
57  return sin(x * 10.) * cos(y * 10.);
58  }
59  static double nablaFunction(const double x, const double y, const double z) {
60  return 200 * sin(x * 10.) * cos(y * 10.);
61  }
62 
63  static int integrationRule(int, int, int p_data) { return 2 * p_data; };
64 
65  MoFEMErrorCode setupProblem();
66  MoFEMErrorCode createCommonData();
67  MoFEMErrorCode bC();
68  MoFEMErrorCode OPs();
69  MoFEMErrorCode kspSolve();
70  MoFEMErrorCode postProcess();
71  MoFEMErrorCode checkResults();
72 
74  struct CommonData {
75  boost::shared_ptr<VectorDouble> approxVals;
78  };
79  boost::shared_ptr<CommonData> commonDataPtr;
80  boost::shared_ptr<std::vector<bool>> boundaryMarker;
81 
82  struct OpError : public DomainEleOp {
83  boost::shared_ptr<CommonData> commonDataPtr;
84  OpError(boost::shared_ptr<CommonData> &common_data_ptr)
85  : DomainEleOp("U", OPROW), commonDataPtr(common_data_ptr) {}
86  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
87  };
88 };
89 
92  CHKERR setupProblem();
93  CHKERR createCommonData();
94  CHKERR bC();
95  CHKERR OPs();
96  CHKERR kspSolve();
97  CHKERR postProcess();
98  CHKERR checkResults();
100 }
101 
102 //! [Set up problem]
105  Simple *simple = mField.getInterface<Simple>();
106  // Add field
107  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
108  CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
109  constexpr int order = 5;
110  CHKERR simple->setFieldOrder("U", order);
111  CHKERR simple->setUp();
113 }
114 //! [Set up problem]
115 
116 //! [Create common data]
119  Simple *simple = mField.getInterface<Simple>();
120  commonDataPtr = boost::make_shared<CommonData>();
121  commonDataPtr->resVec = smartCreateDMVector(simple->getDM());
122  commonDataPtr->L2Vec = createSmartVectorMPI(
123  mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
124  commonDataPtr->approxVals = boost::make_shared<VectorDouble>();
126 }
127 //! [Create common data]
128 
129 //! [Boundary condition]
132 
133  Simple *simple = mField.getInterface<Simple>();
134 
135  auto get_ents_on_mesh_skin = [&]() {
136  Range faces;
137  CHKERR mField.get_moab().get_entities_by_type(0, MBTRI, faces);
138  Skinner skin(&mField.get_moab());
139  Range skin_edges;
140  CHKERR skin.find_skin(0, faces, false, skin_edges);
141  Range skin_verts;
142  CHKERR mField.get_moab().get_connectivity(skin_edges, skin_verts, true);
143  skin_edges.merge(skin_verts);
144  return skin_edges;
145  };
146 
147  auto mark_boundary_dofs = [&](Range &&skin_edges) {
148  auto problem_manager = mField.getInterface<ProblemsManager>();
149  auto marker_ptr = boost::make_shared<std::vector<bool>>();
150  problem_manager->markDofs(simple->getProblemName(), ROW, skin_edges,
151  *marker_ptr);
152  return marker_ptr;
153  };
154 
155  // Get global local vector of marked DOFs. Is global, since is set for all
156  // DOFs on processor. Is local since only DOFs on processor are in the
157  // vector. To access DOFs use local indices.
158  boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
159 
161 }
162 //! [Boundary condition]
163 
164 //! [Push operators to pipeline]
167  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
168  pipeline_mng->getOpDomainLhsPipeline().push_back(
170  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetInvJacH1ForFace(invJac));
171  auto beta = [](const double, const double, const double) { return 1; };
172  pipeline_mng->getOpDomainLhsPipeline().push_back(
173  new OpDomainGradGrad("U", "U", beta, boundaryMarker));
174  pipeline_mng->getOpDomainRhsPipeline().push_back(
175  new OpDomainSource("U", Example::nablaFunction, boundaryMarker));
176  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integrationRule);
177  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integrationRule);
178  pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpBoundaryMass("U", "U", beta));
179  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
180  new OpBoundarySource("U", approxFunction));
181  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integrationRule);
182  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integrationRule);
184 }
185 //! [Push operators to pipeline]
186 
187 //! [Solve]
190  Simple *simple = mField.getInterface<Simple>();
191  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
192  auto solver = pipeline_mng->createKSP();
193  CHKERR KSPSetFromOptions(solver);
194  CHKERR KSPSetUp(solver);
195 
196  auto dm = simple->getDM();
197  auto D = smartCreateDMVector(dm);
198  auto F = smartVectorDuplicate(D);
199 
200  CHKERR KSPSolve(solver, F, D);
201  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
202  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
203  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
205 }
206 
207 //! [Solve]
210  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
211  pipeline_mng->getDomainLhsFE().reset();
212  pipeline_mng->getBoundaryLhsFE().reset();
213  pipeline_mng->getBoundaryRhsFE().reset();
214  auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
215  post_proc_fe->generateReferenceElementMesh();
216  post_proc_fe->addFieldValuesPostProc("U");
217  pipeline_mng->getDomainRhsFE() = post_proc_fe;
218  CHKERR pipeline_mng->loopFiniteElements();
219  CHKERR post_proc_fe->writeFile("out_poisson.h5m");
221 }
222 //! [Postprocess results]
223 
224 //! [Solve]
227  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
228  pipeline_mng->getDomainLhsFE().reset();
229  pipeline_mng->getDomainRhsFE().reset();
230  pipeline_mng->getBoundaryLhsFE().reset();
231  pipeline_mng->getBoundaryRhsFE().reset();
232  pipeline_mng->getOpDomainRhsPipeline().clear();
233  pipeline_mng->getOpDomainRhsPipeline().push_back(
234  new OpCalculateScalarFieldValues("U", commonDataPtr->approxVals));
235  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpError(commonDataPtr));
236  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integrationRule);
237  CHKERR pipeline_mng->loopFiniteElements();
238  CHKERR VecAssemblyBegin(commonDataPtr->L2Vec);
239  CHKERR VecAssemblyEnd(commonDataPtr->L2Vec);
240  CHKERR VecAssemblyBegin(commonDataPtr->resVec);
241  CHKERR VecAssemblyEnd(commonDataPtr->resVec);
242  double nrm2;
243  CHKERR VecNorm(commonDataPtr->resVec, NORM_2, &nrm2);
244  const double *array;
245  CHKERR VecGetArrayRead(commonDataPtr->L2Vec, &array);
246  if (mField.get_comm_rank() == 0)
247  PetscPrintf(PETSC_COMM_SELF, "Error %6.4e Vec norm %6.4e\n", sqrt(array[0]),
248  nrm2);
249  CHKERR VecRestoreArrayRead(commonDataPtr->L2Vec, &array);
250  constexpr double eps = 1e-8;
251  if (nrm2 > eps)
252  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
253  "Not converged solution");
255 }
256 //! [Solver]
257 
258 int main(int argc, char *argv[]) {
259 
260  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
261 
262  try {
263 
264  //! [Register MoFEM discrete manager in PETSc]
265  DMType dm_name = "DMMOFEM";
266  CHKERR DMRegister_MoFEM(dm_name);
267  //! [Register MoFEM discrete manager in PETSc
268 
269  //! [Create MoAB]
270  moab::Core mb_instance; ///< mesh database
271  moab::Interface &moab = mb_instance; ///< mesh database interface
272  //! [Create MoAB]
273 
274  //! [Create MoFEM]
275  MoFEM::Core core(moab); ///< finite element database
276  MoFEM::Interface &m_field = core; ///< finite element database insterface
277  //! [Create MoFEM]
278 
279  //! [Load mesh]
280  Simple *simple = m_field.getInterface<Simple>();
281  CHKERR simple->getOptions();
282  CHKERR simple->loadFile("");
283  //! [Load mesh]
284 
285  //! [Example]
286  Example ex(m_field);
287  CHKERR ex.runProblem();
288  //! [Example]
289  }
290  CATCH_ERRORS;
291 
293 }
294 
296  EntData &data) {
298 
299  if (const size_t nb_dofs = data.getIndices().size()) {
300 
301  const int nb_integration_pts = getGaussPts().size2();
302  auto t_w = getFTensor0IntegrationWeight();
303  auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
304  auto t_coords = getFTensor1CoordsAtGaussPts();
305 
306  VectorDouble nf(nb_dofs, false);
307  nf.clear();
308 
310  const double volume = getMeasure();
311 
312  auto t_row_base = data.getFTensor0N();
313  double error = 0;
314  for (int gg = 0; gg != nb_integration_pts; ++gg) {
315 
316  const double alpha = t_w * volume;
317  double diff = t_val - Example::approxFunction(t_coords(0), t_coords(1),
318  t_coords(2));
319  error += alpha * pow(diff, 2);
320 
321  for (size_t r = 0; r != nb_dofs; ++r) {
322  nf[r] += alpha * t_row_base * diff;
323  ++t_row_base;
324  }
325 
326  ++t_w;
327  ++t_val;
328  ++t_coords;
329  }
330 
331  const int index = 0;
332  CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
333  CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
334  }
335 
337 }
boost::shared_ptr< CommonData > commonDataPtr
SmartPetscObj< Vec > resVec
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:143
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 setupProblem()
[Set up problem]
MoFEMErrorCode kspSolve()
[Push operators to pipeline]
CoreTmp< 0 > Core
Definition: Core.hpp:1129
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
MoFEM::Interface & mField
const double D
diffusivity
[Example]
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
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
boost::shared_ptr< std::vector< bool > > boundaryMarker
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:485
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
Example(MoFEM::Interface &m_field)
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
SmartPetscObj< Vec > L2Vec
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.
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()
[Run problem]
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.
constexpr int order
MoFEMErrorCode checkResults()
[Postprocess results]
Calculate inverse of jacobian for face element.
MoFEMErrorCode bC()
[Create common data]
static Index< 'i', 3 > i
[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 int integrationRule(int, int, int p_data)
PipelineManager interface.
#define CHKERR
Inline error check.
Definition: definitions.h:604
boost::shared_ptr< FEMethod > & getBoundaryLhsFE()
const double r
rate factor
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
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< VectorDouble > approxVals
MoFEMErrorCode postProcess()
[Solve]
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:1943
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.
Core (interface) class.
Definition: Core.hpp:77
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
continuous field
Definition: definitions.h:177
static const double eps
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
Get value at integration points for scalar field.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
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.
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...