v0.9.2
lesson6_radiation.cpp
Go to the documentation of this file.
1 /**
2  * \file lesson6_radiation.cpp
3  * \example lesson6_radiation.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 
43 
44 constexpr double emissivity = 1;
45 constexpr double boltzmann_constant = 5.670367e-2;
46 constexpr double Beta = emissivity * boltzmann_constant;
47 
48 constexpr double T_ambient = 3;
49 constexpr double Flux = -1.0e6;
50 
51 struct Example {
52 
53  Example(MoFEM::Interface &m_field) : mField(m_field) {}
54 
55  MoFEMErrorCode runProblem();
56 
57 private:
58  MoFEM::Interface &mField;
59 
60  static int integrationRule(int, int, int p_data) { return 2 * p_data; };
61 
62  MoFEMErrorCode setUP();
63  MoFEMErrorCode createCommonData();
64  MoFEMErrorCode bC();
65  MoFEMErrorCode OPs();
66  MoFEMErrorCode kspSolve();
67  MoFEMErrorCode postProcess();
68  MoFEMErrorCode checkResults();
69 
71  boost::shared_ptr<VectorDouble> approxVals;
72  boost::shared_ptr<MatrixDouble> approxGradVals;
73 
74  struct OpRadiationLhs : public OpFaceBase {
75 
76  private:
77  boost::shared_ptr<VectorDouble> approxVals;
78  FTensor::Index<'i', 2> i; ///< summit Index
79 
80  public:
81  OpRadiationLhs(boost::shared_ptr<VectorDouble> &approx_vals)
82  : OpFaceBase("U", "U", OpFaceBase::OPROWCOL), approxVals(approx_vals) {}
83 
84  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
85  };
86 
87  struct OpRadiationRhs : public OpFaceBase {
88 
89  private:
90  boost::shared_ptr<VectorDouble> approxVals;
91  FTensor::Index<'i', 2> i; ///< summit Index
92 
93  public:
94  OpRadiationRhs(boost::shared_ptr<VectorDouble> &approx_vals)
95  : OpFaceBase("U", "U", OpFaceBase::OPROW), approxVals(approx_vals) {}
96 
97  MoFEMErrorCode iNtegrate(EntData &row_data);
98  };
99 
100  struct OpFluxRhs : public OpFaceBase {
101 
102  private:
103  FTensor::Index<'i', 2> i; ///< summit Index
104 
105  public:
106  OpFluxRhs() : OpFaceBase("U", "U", OpFaceBase::OPROW) {}
107 
108  MoFEMErrorCode iNtegrate(EntData &row_data);
109  };
110 };
111 
114  CHKERR setUP();
115  CHKERR createCommonData();
116  CHKERR bC();
117  CHKERR OPs();
118  CHKERR kspSolve();
119  CHKERR postProcess();
120  CHKERR checkResults();
122 }
123 
124 //! [Set up problem]
127  Simple *simple = mField.getInterface<Simple>();
128  // Add field
129  CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
130  CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
131  constexpr int order = 3;
132  CHKERR simple->setFieldOrder("U", order);
133  CHKERR simple->setUp();
135 }
136 //! [Set up problem]
137 
138 //! [Create common data]
141  approxVals = boost::make_shared<VectorDouble>();
142  approxGradVals = boost::make_shared<MatrixDouble>();
144 }
145 //! [Create common data]
146 
147 //! [Boundary condition]
150  // Set initial values
151  auto set_initial_temperature = [&](VectorAdaptor &&field_data, double *xcoord,
152  double *ycoord, double *zcoord) {
154  field_data[0] = T_ambient;
156  };
157  FieldBlas *field_blas;
158  CHKERR mField.getInterface(field_blas);
159  CHKERR field_blas->setVertexDofs(set_initial_temperature, "U");
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 r, const double, const double) {
172  return 2e3 * 2 * M_PI * r;
173  };
174  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpDomainGradGrad("U", "U", beta));
175  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integrationRule);
176 
177  pipeline_mng->getOpDomainRhsPipeline().push_back(
179  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSetInvJacH1ForFace(invJac));
180  pipeline_mng->getOpDomainRhsPipeline().push_back(
181  new OpCalculateScalarFieldGradient<2>("U", approxGradVals));
182  pipeline_mng->getOpDomainRhsPipeline().push_back(
183  new OpVolGradGradResidual("U", beta, approxGradVals));
184  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integrationRule);
185 
186  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
187  new OpCalculateScalarFieldValues("U", approxVals));
188  pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpRadiationRhs(approxVals));
189  pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpFluxRhs());
190  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integrationRule);
191 
192  pipeline_mng->getOpBoundaryLhsPipeline().push_back(
193  new OpCalculateScalarFieldValues("U", approxVals));
194  pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpRadiationLhs(approxVals));
195  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integrationRule);
197 }
198 //! [Push operators to pipeline]
199 
200 //! [Solve]
203  Simple *simple = mField.getInterface<Simple>();
204  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
205  auto ts = pipeline_mng->createTS();
206 
207  double ftime = 1;
208  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
209  CHKERR TSSetFromOptions(ts);
210  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
211 
212  auto T = smartCreateDMVector(simple->getDM());
213  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
214  SCATTER_FORWARD);
215 
216  CHKERR TSSolve(ts, T);
217  CHKERR TSGetTime(ts, &ftime);
218 
219  PetscInt steps, snesfails, rejects, nonlinits, linits;
220  CHKERR TSGetTimeStepNumber(ts, &steps);
221  CHKERR TSGetSNESFailures(ts, &snesfails);
222  CHKERR TSGetStepRejections(ts, &rejects);
223  CHKERR TSGetSNESIterations(ts, &nonlinits);
224  CHKERR TSGetKSPIterations(ts, &linits);
225  PetscPrintf(PETSC_COMM_WORLD,
226  "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
227  "%D, linits %D\n",
228  steps, rejects, snesfails, ftime, nonlinits, linits);
229 
231 }
232 //! [Solve]
233 
234 //! [Postprocess results]
237  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
238  pipeline_mng->getDomainLhsFE().reset();
239  pipeline_mng->getBoundaryLhsFE().reset();
240  pipeline_mng->getBoundaryRhsFE().reset();
241  auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
242  post_proc_fe->generateReferenceElementMesh();
243  post_proc_fe->addFieldValuesPostProc("U");
244  pipeline_mng->getDomainRhsFE() = post_proc_fe;
245  CHKERR pipeline_mng->loopFiniteElements();
246  CHKERR post_proc_fe->writeFile("out_radiation.h5m");
248 }
249 //! [Postprocess results]
250 
251 //! [Check results]
253 //! [Check results]
254 
255 int main(int argc, char *argv[]) {
256 
257  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
258 
259  try {
260 
261  //! [Register MoFEM discrete manager in PETSc]
262  DMType dm_name = "DMMOFEM";
263  CHKERR DMRegister_MoFEM(dm_name);
264  //! [Register MoFEM discrete manager in PETSc
265 
266  //! [Create MoAB]
267  moab::Core mb_instance; ///< mesh database
268  moab::Interface &moab = mb_instance; ///< mesh database interface
269  //! [Create MoAB]
270 
271  //! [Create MoFEM]
272  MoFEM::Core core(moab); ///< finite element database
273  MoFEM::Interface &m_field = core; ///< finite element database insterface
274  //! [Create MoFEM]
275 
276  //! [Load mesh]
277  Simple *simple = m_field.getInterface<Simple>();
278  CHKERR simple->getOptions();
279  CHKERR simple->loadFile("");
280  //! [Load mesh]
281 
282  //! [Example]
283  Example ex(m_field);
284  CHKERR ex.runProblem();
285  //! [Example]
286  }
287  CATCH_ERRORS;
288 
290 }
291 
292 //! [Radiation Lhs]
294  EntData &col_data) {
296  // get element volume
297  const double vol = OpFaceBase::getMeasure();
298  // get integration weights
299  auto t_w = OpFaceBase::getFTensor0IntegrationWeight();
300  // get base function gradient on rows
301  auto t_row_base = row_data.getFTensor0N();
302  // gat temperature at integration points
303  auto t_val = getFTensor0FromVec(*(approxVals));
304  // get coordinate at integration points
305  auto t_coords = OpFaceBase::getFTensor1CoordsAtGaussPts();
306 
307  // loop over integration points
308  for (int gg = 0; gg != OpFaceBase::nbIntegrationPts; gg++) {
309 
310  // Cylinder radius
311  const double r_cylinder = t_coords(0);
312 
313  // take into account Jacobean
314  const double alpha = t_w * vol * Beta * (2 * M_PI * r_cylinder);
315  // loop over rows base functions
316  for (int rr = 0; rr != OpFaceBase::nbRows; ++rr) {
317  auto t_col_base = col_data.getFTensor0N(gg, 0);
318  // loop over columns
319  for (int cc = 0; cc != OpFaceBase::nbCols; cc++) {
320  if (std::abs(t_coords(0)) > std::numeric_limits<double>::epsilon()) {
321  locMat(rr, cc) += alpha * t_row_base * t_col_base * 4 * pow(t_val, 3);
322  }
323  ++t_col_base;
324  }
325  ++t_row_base;
326  }
327 
328  ++t_val;
329  ++t_coords;
330  ++t_w; // move to another integration weight
331  }
333 }
334 //! [Radiation Lhs]
335 
336 //! [Radiation Lhs]
339  // get element volume
340  const double vol = OpFaceBase::getMeasure();
341  // get integration weights
342  auto t_w = OpFaceBase::getFTensor0IntegrationWeight();
343  // get base function gradient on rows
344  auto t_row_base = row_data.getFTensor0N();
345  // gat temperature at integration points
346  auto t_val = getFTensor0FromVec(*(approxVals));
347  // get coordinate at integration points
348  auto t_coords = OpFaceBase::getFTensor1CoordsAtGaussPts();
349 
350  // loop over integration points
351  for (int gg = 0; gg != OpFaceBase::nbIntegrationPts; gg++) {
352 
353  // Cylinder radius
354  const double r_cylinder = t_coords(0);
355 
356  // take into account Jacobean
357  const double alpha = t_w * vol * Beta * (2 * M_PI * r_cylinder);
358  // loop over rows base functions
359  for (int rr = 0; rr != OpFaceBase::nbRows; ++rr) {
360  if (std::abs(t_coords(0)) > std::numeric_limits<double>::epsilon()) {
361  OpFaceBase::locF[rr] +=
362  alpha * t_row_base * (pow(t_val, 4) - pow(T_ambient, 4));
363  }
364  ++t_row_base;
365  }
366  ++t_coords;
367  ++t_val;
368  ++t_w; // move to another integration weight
369  }
370 
372 }
373 //! [Radiation Lhs]
374 
375 //! [Flux Rhs]
378  // get element volume
379  const double vol = OpFaceBase::getMeasure();
380  // get integration weights
381  auto t_w = OpFaceBase::getFTensor0IntegrationWeight();
382  // get base function gradient on rows
383  auto t_row_base = row_data.getFTensor0N();
384  // get coordinate at integration points
385  auto t_coords = OpFaceBase::getFTensor1CoordsAtGaussPts();
386  // get time
387  const double time = OpFaceBase::getFEMethod()->ts_t;
388 
389  // loop over integration points
390  for (int gg = 0; gg != OpFaceBase::nbIntegrationPts; gg++) {
391 
392  // Cylinder radius
393  const double r_cylinder = t_coords(0);
394 
395  const double r = sqrt(t_coords(i) * t_coords(i));
396  const double c = std::abs(t_coords(0)) / r;
397 
398  // take into account Jacobean
399  const double alpha = t_w * vol * (2 * M_PI * r_cylinder);
400  // loop over rows base functions
401  for (int rr = 0; rr != OpFaceBase::nbRows; ++rr) {
402  OpFaceBase::locF[rr] += alpha * t_row_base * c * Flux * time;
403  ++t_row_base;
404  }
405  ++t_coords;
406  ++t_w; // move to another integration weight
407  }
408 
410 }
411 //! [Flux Rhs]
const double T
FTensor::Index< 'i', 2 > i
summit Index
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.
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Integrate B^T D B operator.
Deprecated interface functions.
OpRadiationLhs(boost::shared_ptr< VectorDouble > &approx_vals)
MoFEMErrorCode createCommonData()
[Set up problem]
constexpr double emissivity
MoFEMErrorCode kspSolve()
[Push operators to pipeline]
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
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
SmartPetscObj< TS > createTS(SmartPetscObj< DM > dm=nullptr)
Create TS (time) solver.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
FTensor::Index< 'i', 2 > i
summit Index
MatrixDouble invJac
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
int nbCols
number if dof on column
Definition: BaseOps.hpp:57
Example(MoFEM::Interface &m_field)
Core (interface) class.
Definition: Core.hpp:50
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
[Radiation Lhs]
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
boost::shared_ptr< VectorDouble > approxVals
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
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...
constexpr double Flux
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:939
MoFEMErrorCode iNtegrate(EntData &row_data)
[Radiation Lhs]
constexpr double boltzmann_constant
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]
MoFEMErrorCode checkResults()
[Print results]
Calculate inverse of jacobian for face element.
MoFEMErrorCode bC()
[Create common data]
VectorDouble locF
local entity vector
Definition: BaseOps.hpp:61
Basic algebra on fields.
Definition: FieldBlas.hpp:34
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
DataForcesAndSourcesCore::EntData EntData
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:108
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
boost::shared_ptr< VectorDouble > approxVals
static int integrationRule(int, int, int p_data)
constexpr double Beta
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.
boost::shared_ptr< VectorDouble > approxVals
Transform local reference derivatives of shape functions to global derivatives.
constexpr int order
constexpr double T_ambient
int nbIntegrationPts
number of integration points
Definition: BaseOps.hpp:58
boost::shared_ptr< MatrixDouble > approxGradVals
MoFEMErrorCode postProcess()
[Integrate]
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
boost::ptr_vector< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
[Source operator]
Definition: BaseOps.hpp:137
static char help[]
int nbRows
number of dofs on rows
Definition: BaseOps.hpp:56
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
#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
OpRadiationRhs(boost::shared_ptr< VectorDouble > &approx_vals)
continuous field
Definition: definitions.h:177
FTensor::Index< 'i', 2 > i
summit Index
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode iNtegrate(EntData &row_data)
[Radiation Lhs]
OpTools< DomainEleOp >::OpGradGradResidual< 2 > OpVolGradGradResidual
Get value at integration points for scalar field.
OpTools< DomainEleOp >::OpGradGrad< 2 > OpDomainGradGrad
int main(int argc, char *argv[])
[Check results]
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
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
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user functionExample:
Definition: FieldBlas.cpp:172
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