v0.13.1
Loading...
Searching...
No Matches
PoissonOperators.hpp
Go to the documentation of this file.
1/**
2 * \file PoissonOperators.hpp
3 * \example PoissonOperators.hpp
4 *
5 */
6
7#ifndef __POISSONOPERATORS_HPP__
8#define __POISSONOPERATORS_HPP__
9
10namespace PoissonExample {
11
12using PostProcFE = PostProcBrokenMeshInMoab<VolumeElementForcesAndSourcesCore>;
13
14/**
15 * \brief Calculate the grad-grad operator and assemble matrix
16 *
17 * Calculate
18 * \f[
19 * \mathbf{K}=\int_\Omega \nabla \boldsymbol\phi \cdot \nabla \boldsymbol\phi
20 * \textrm{d}\Omega \f] and assemble to global matrix.
21 *
22 * This operator is executed on element for each unique combination of entities.
23 *
24 */
25struct OpK : public VolumeElementForcesAndSourcesCore::UserDataOperator {
26
27 OpK(bool symm = true)
28 : VolumeElementForcesAndSourcesCore::UserDataOperator("U", "U", OPROWCOL,
29 symm) {}
30
31 /**
32 * \brief Do calculations for give operator
33 * @param row_side row side number (local number) of entity on element
34 * @param col_side column side number (local number) of entity on element
35 * @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
36 * @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
37 * @param row_data data for row
38 * @param col_data data for column
39 * @return error code
40 */
41 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
42 EntityType col_type,
43 EntitiesFieldData::EntData &row_data,
44 EntitiesFieldData::EntData &col_data) {
46 // get number of dofs on row
47 nbRows = row_data.getIndices().size();
48 // if no dofs on row, exit that work, nothing to do here
49 if (!nbRows)
51 // get number of dofs on column
52 nbCols = col_data.getIndices().size();
53 // if no dofs on Columbia, exit nothing to do here
54 if (!nbCols)
56 // get number of integration points
57 nbIntegrationPts = getGaussPts().size2();
58 // check if entity block is on matrix diagonal
59 if (row_side == col_side && row_type == col_type) {
60 isDiag = true; // yes, it is
61 } else {
62 isDiag = false;
63 }
64 // integrate local matrix for entity block
65 CHKERR iNtegrate(row_data, col_data);
66 // assemble local matrix
67 CHKERR aSsemble(row_data, col_data);
69 }
70
71protected:
72 ///< error code
73
74 int nbRows; ///< number of dofs on rows
75 int nbCols; ///< number if dof on column
76 int nbIntegrationPts; ///< number of integration points
77 bool isDiag; ///< true if this block is on diagonal
78
79 FTensor::Index<'i', 3> i; ///< summit Index
80 MatrixDouble locMat; ///< local entity block matrix
81
82 /**
83 * \brief Integrate grad-grad operator
84 * @param row_data row data (consist base functions on row entity)
85 * @param col_data column data (consist base functions on column entity)
86 * @return error code
87 */
88 virtual MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
89 EntitiesFieldData::EntData &col_data) {
91 // set size of local entity bock
92 locMat.resize(nbRows, nbCols, false);
93 // clear matrix
94 locMat.clear();
95 // get element volume
96 double vol = getVolume();
97 // get integration weights
98 auto t_w = getFTensor0IntegrationWeight();
99 // get base function gradient on rows
100 auto t_row_grad = row_data.getFTensor1DiffN<3>();
101 // loop over integration points
102 for (int gg = 0; gg != nbIntegrationPts; gg++) {
103 // take into account Jacobean
104 const double alpha = t_w * vol;
105 // noalias(locMat) +=
106 // alpha*prod(row_data.getDiffN(gg),trans(col_data.getDiffN(gg))); take
107 // fist element to local matrix
108 FTensor::Tensor0<double *> a(&*locMat.data().begin());
109 // loop over rows base functions
110 for (int rr = 0; rr != nbRows; rr++) {
111 // get column base functions gradient at gauss point gg
112 auto t_col_grad = col_data.getFTensor1DiffN<3>(gg, 0);
113 // loop over columns
114 for (int cc = 0; cc != nbCols; cc++) {
115 // calculate element of local matrix
116 a += alpha * (t_row_grad(i) * t_col_grad(i));
117 ++t_col_grad; // move to another gradient of base function on column
118 ++a; // move to another element of local matrix in column
119 }
120 ++t_row_grad; // move to another element of gradient of base function on
121 // row
122 }
123 ++t_w; // move to another integration weight
124 }
126 }
127
128 /**
129 * \brief Assemble local entity block matrix
130 * @param row_data row data (consist base functions on row entity)
131 * @param col_data column data (consist base functions on column entity)
132 * @return error code
133 */
134 virtual MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data,
135 EntitiesFieldData::EntData &col_data) {
137 // get pointer to first global index on row
138 const int *row_indices = &*row_data.getIndices().data().begin();
139 // get pointer to first global index on column
140 const int *col_indices = &*col_data.getIndices().data().begin();
141 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
142 : getFEMethod()->snes_B;
143 // assemble local matrix
144 CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
145 &*locMat.data().begin(), ADD_VALUES);
146
147 if (!isDiag && sYmm) {
148 // if not diagonal term and since global matrix is symmetric assemble
149 // transpose term.
150 locMat = trans(locMat);
151 CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
152 &*locMat.data().begin(), ADD_VALUES);
153 }
155 }
156};
157
158/**
159 * \brief template class for integration oh the right hand side
160 */
161template <typename OPBASE> struct OpBaseRhs : public OPBASE {
162
163 OpBaseRhs(const std::string field_name) : OPBASE(field_name, OPBASE::OPROW) {}
164
165 /**
166 * \brief This function is called by finite element
167 *
168 * Do work is composed from two operations, integrate and assembly. Also,
169 * it set values nbRows, and nbIntegrationPts.
170 *
171 */
172 MoFEMErrorCode doWork(int row_side, EntityType row_type,
173 EntitiesFieldData::EntData &row_data) {
175 // get number of dofs on row
176 nbRows = row_data.getIndices().size();
177 if (!nbRows)
179 // get number of integration points
180 nbIntegrationPts = OPBASE::getGaussPts().size2();
181 // integrate local vector
182 CHKERR iNtegrate(row_data);
183 // assemble local vector
184 CHKERR aSsemble(row_data);
186 }
187
188 /**
189 * \brief Class dedicated to integrate operator
190 * @param data entity data on element row
191 * @return error code
192 */
193 virtual MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) = 0;
194
195 /**
196 * \brief Class dedicated to assemble operator to global system vector
197 * @param data entity data (indices, base functions, etc. ) on element row
198 * @return error code
199 */
200 virtual MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data) = 0;
201
202protected:
203 ///< error code
204 int nbRows; ///< number of dofs on row
205 int nbIntegrationPts; ///< number of integration points
206};
207
208/**
209 * \brief Operator calculate source term,
210 *
211 * \f[
212 * \mathbf{F} = \int_\Omega \boldsymbol\phi f \textrm{d}\Omega
213 * \f]
214 *
215 */
216struct OpF
217 : public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
218
219 typedef boost::function<double(const double, const double, const double)>
221
222 OpF(FSource f_source)
223 : OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator>("U"),
224 fSource(f_source) {}
225
226protected:
231
232 VectorDouble locVec;
233
234 /**
235 * \brief Integrate local entity vector
236 * @param data entity data on element row
237 * @return error code
238 */
239 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
241 // set size of local vector
242 locVec.resize(nbRows, false);
243 // clear local entity vector
244 locVec.clear();
245 // get finite element volume
246 double vol = getVolume();
247 // get integration weights
248 auto t_w = getFTensor0IntegrationWeight();
249 // get base functions on entity
250 auto t_v = data.getFTensor0N();
251 // get coordinates at integration points
252 auto t_coords = getFTensor1CoordsAtGaussPts();
253 // loop over all integration points
254 for (int gg = 0; gg != nbIntegrationPts; gg++) {
255 // evaluate constant term
256 const double alpha =
257 vol * t_w * fSource(t_coords(NX), t_coords(NY), t_coords(NZ));
258 // get element of local vector
260 &*locVec.data().begin());
261 // loop over base functions
262 for (int rr = 0; rr != nbRows; rr++) {
263 // add to local vector source term
264 t_a -= alpha * t_v;
265 ++t_a; // move to next element of local vector
266 ++t_v; // move to next base function
267 }
268 ++t_w; // move to next integration weight
269 ++t_coords; // move to next physical coordinates at integration point
270 }
272 }
273
274 /**
275 * \brief assemble local entity vector to the global right hand side
276 * @param data entity data, i.e. global indices of local vector
277 * @return error code
278 */
279 MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data) {
281 // get global indices of local vector
282 const int *indices = &*data.getIndices().data().begin();
283 // get values from local vector
284 const double *vals = &*locVec.data().begin();
285 Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
286 : getFEMethod()->snes_f;
287 // assemble vector
288 CHKERR VecSetValues(f, nbRows, indices, vals, ADD_VALUES);
290 }
291};
292
293/**
294 * \brief Calculate constrains matrix
295 *
296 * \f[
297 * \mathbf{C} = \int_{\partial\Omega} \boldsymbol\psi \boldsymbol\phi
298 * \textrm{d}\partial\Omega \f] where \f$\lambda \f$ is base function on
299 * boundary
300 *
301 */
302struct OpC : public FaceElementForcesAndSourcesCore::UserDataOperator {
303
304 OpC(const bool assemble_transpose)
306 false),
307 assembleTranspose(assemble_transpose) {}
308
309 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
310 EntityType col_type,
311 EntitiesFieldData::EntData &row_data,
312 EntitiesFieldData::EntData &col_data) {
314 // get number of dofs on row
315 nbRows = row_data.getIndices().size();
316 // exit here if no dofs on row, nothing to do
317 if (!nbRows)
319 // get number of dofs on column,
320 nbCols = col_data.getIndices().size();
321 // exit here if no dofs on roe, nothing to do
322 if (!nbCols)
324 // get number of integration points
325 nbIntegrationPts = getGaussPts().size2();
326 // integrate local constrains matrix
327 CHKERR iNtegrate(row_data, col_data);
328 // assemble local constrains matrix
329 CHKERR aSsemble(row_data, col_data);
331 }
332
333private:
334 ///< error code
335
336 int nbRows; ///< number of dofs on row
337 int nbCols; ///< number of dofs on column
338 int nbIntegrationPts; ///< number of integration points
339 const bool assembleTranspose; ///< assemble transpose, i.e. CT if set to true
340
341 MatrixDouble locMat; ///< local constrains matrix
342
343 /** \brief Integrate local constrains matrix
344 */
345 inline MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
346 EntitiesFieldData::EntData &col_data) {
348 // set size of local constrains matrix
349 locMat.resize(nbRows, nbCols, false);
350 // clear matrix
351 locMat.clear();
352 // get area of element
353 const double area = getArea();
354 // get integration weights
355 auto t_w = getFTensor0IntegrationWeight();
356 // get base functions on entity
357 auto t_row = row_data.getFTensor0N();
358 // run over integration points
359 for (int gg = 0; gg != nbIntegrationPts; gg++) {
360 const double alpha = area * t_w;
361 // get element of local matrix
363 &*locMat.data().begin());
364 // run over base functions on rows
365 for (int rr = 0; rr != nbRows; rr++) {
366 // get first base functions on column for integration point gg
367 auto t_col = col_data.getFTensor0N(gg, 0);
368 // run over base function on column
369 for (int cc = 0; cc != nbCols; cc++) {
370 // integrate element of constrains matrix
371 c += alpha * t_row * t_col;
372 ++t_col; // move to next base function on column
373 ++c; // move to next element of local matrix
374 }
375 ++t_row; // move to next base function on row
376 }
377 ++t_w; // move to next integrate weight
378 }
380 }
381
382 /**
383 * \brief integrate local constrains matrix
384 */
385 inline MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data,
386 EntitiesFieldData::EntData &col_data) {
388 // get indices on row
389 const int *row_indices = &*row_data.getIndices().data().begin();
390 // get indices on column
391 const int *col_indices = &*col_data.getIndices().data().begin();
392 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
393 : getFEMethod()->snes_B;
394 // assemble local matrix
395 CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
396 &*locMat.data().begin(), ADD_VALUES);
397 // cerr << locMat << endl;
398 if (assembleTranspose) {
399 // assemble transpose of local matrix
400 locMat = trans(locMat);
401 CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
402 &*locMat.data().begin(), ADD_VALUES);
403 }
405 }
406};
407
408/**
409 * \brief Assemble constrains vector
410 *
411 * \f[
412 * \mathbf{g} = \int_{\partial\Omega} \boldsymbol\psi \overline{u}
413 * \textrm{d}\partial\Omega \f]
414 *
415 */
416struct Op_g
417 : public OpBaseRhs<FaceElementForcesAndSourcesCore::UserDataOperator> {
418
419 typedef boost::function<double(const double, const double, const double)>
421
422 Op_g(FVal f_value, const string field_name = "L", const double beta = 1)
424 field_name),
425 fValue(f_value), bEta(beta) {}
426
427protected:
428 FTensor::Number<0> NX; ///< x-direction index
429 FTensor::Number<1> NY; ///< y-direction index
430 FTensor::Number<2> NZ; ///< z-direction index
431 FVal fValue; ///< Function pointer evaluating values of "U" at the boundary
432
433 VectorDouble locVec;
434 const double bEta;
435
436 /**
437 * \brief Integrate local constrains vector
438 */
439 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
441 // set size to local vector
442 locVec.resize(nbRows, false);
443 // clear local vector
444 locVec.clear();
445 // get face area
446 const double area = getArea() * bEta;
447 // get integration weight
448 auto t_w = getFTensor0IntegrationWeight();
449 // get base function
450 auto t_l = data.getFTensor0N();
451 // get coordinates at integration point
452 auto t_coords = getFTensor1CoordsAtGaussPts();
453 // make loop over integration points
454 for (int gg = 0; gg != nbIntegrationPts; gg++) {
455 // evaluate function on boundary and scale it by area and integration
456 // weight
457 double alpha =
458 area * t_w * fValue(t_coords(NX), t_coords(NY), t_coords(NZ));
459 // get element of vector
461 &*locVec.data().begin());
462 //
463 for (int rr = 0; rr != nbRows; rr++) {
464 t_a += alpha * t_l;
465 ++t_a;
466 ++t_l;
467 }
468 ++t_w;
469 ++t_coords;
470 }
472 }
473
474 /**
475 * \brief assemble constrains vectors
476 */
477 MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data) {
479 const int *indices = &*data.getIndices().data().begin();
480 const double *vals = &*locVec.data().begin();
481 Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
482 : getFEMethod()->snes_f;
483 CHKERR VecSetValues(f, nbRows, indices, &*vals, ADD_VALUES);
485 }
486};
487
488/**
489 * \brief Evaluate error
490 */
492 : public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
493
494 typedef boost::function<double(const double, const double, const double)>
496 typedef boost::function<FTensor::Tensor1<double, 3>(
497 const double, const double, const double)>
499
500 OpError(UVal u_value, GVal g_value,
501 boost::shared_ptr<VectorDouble> &field_vals,
502 boost::shared_ptr<MatrixDouble> &grad_vals, Vec global_error)
503 : OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator>("ERROR"),
504 globalError(global_error), uValue(u_value), gValue(g_value),
505 fieldVals(field_vals), gradVals(grad_vals) {}
506
507 MoFEMErrorCode doWork(int row_side, EntityType row_type,
508 EntitiesFieldData::EntData &row_data) {
510 nbRows = row_data.getFieldData().size();
511 if (!nbRows)
513 nbIntegrationPts = getGaussPts().size2();
514 CHKERR iNtegrate(row_data);
515 CHKERR aSsemble(row_data);
517 }
518
519private:
520 Vec globalError; ///< ghost vector with global (integrated over volume) error
521
526 UVal uValue; ///< function with exact solution
527 GVal gValue; ///< function with exact solution for gradient
528
529 boost::shared_ptr<VectorDouble> fieldVals;
530 boost::shared_ptr<MatrixDouble> gradVals;
531
532 /**
533 * \brief Integrate error
534 */
535 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
537 // clear field dofs
538 data.getFieldData().clear();
539 // get volume of element
540 const double vol = getVolume();
541 // get integration weight
542 auto t_w = getFTensor0IntegrationWeight();
543 // get solution at integration point
544 auto t_u = getFTensor0FromVec(*fieldVals);
545 // get solution at integration point
546 auto t_grad = getFTensor1FromMat<3>(*gradVals);
547 // get coordinates at integration point
548 auto t_coords = getFTensor1CoordsAtGaussPts();
549 // keep exact gradient and error or gradient
550 FTensor::Tensor1<double, 3> t_exact_grad, t_error_grad;
551 // integrate over
552 for (int gg = 0; gg != nbIntegrationPts; gg++) {
553 double alpha = vol * t_w;
554 // evaluate exact value
555 double exact_u = uValue(t_coords(NX), t_coords(NY), t_coords(NZ));
556 // evaluate exact gradient
557 t_exact_grad = gValue(t_coords(NX), t_coords(NY), t_coords(NZ));
558 // calculate gradient errro
559 t_error_grad(i) = t_grad(i) - t_exact_grad(i);
560 // error
561 double error = pow(t_u - exact_u, 2) + t_error_grad(i) * t_error_grad(i);
562 // iterate over base functions
563 data.getFieldData()[0] += alpha * error;
564 ++t_w; // move to next integration point
565 ++t_u; // next value of function at integration point
566 ++t_grad; // next gradient at integration point
567 ++t_coords; // next coordinate at integration point
568 }
570 }
571
572 /**
573 * \brief Assemble error
574 */
575 MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data) {
577 // set error on mesh
578 data.getFieldDofs()[0]->getFieldData() = sqrt(data.getFieldData()[0]);
579 // assemble vector to global error
580 CHKERR VecSetValue(globalError, 0, data.getFieldData()[0], ADD_VALUES);
582 }
583};
584
585struct OpKt : public OpK {
586
587 OpKt(boost::function<double(const double)> a,
588 boost::function<double(const double)> diff_a,
589 boost::shared_ptr<VectorDouble> &field_vals,
590 boost::shared_ptr<MatrixDouble> &grad_vals)
591 : OpK(false), A(a), diffA(diff_a), fieldVals(field_vals),
592 gradVals(grad_vals) {}
593
594protected:
595 /**
596 * \brief Integrate grad-grad operator
597 * @param row_data row data (consist base functions on row entity)
598 * @param col_data column data (consist base functions on column entity)
599 * @return error code
600 */
601 inline MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
602 EntitiesFieldData::EntData &col_data) {
604 // set size of local entity bock
605 locMat.resize(nbRows, nbCols, false);
606 // clear matrix
607 locMat.clear();
608 // get element volume
609 double vol = getVolume();
610 // get integration weights
611 auto t_w = getFTensor0IntegrationWeight();
612 // get solution at integration point
613 auto t_u = getFTensor0FromVec(*fieldVals);
614 // get solution at integration point
615 auto t_grad = getFTensor1FromMat<3>(*gradVals);
616 // get base function gradient on rows
617 auto t_row_grad = row_data.getFTensor1DiffN<3>();
618 // loop over integration points
619 for (int gg = 0; gg != nbIntegrationPts; gg++) {
620 // take into account Jacobian
621 const double alpha = t_w * vol;
622 const double beta = alpha * A(t_u);
624 t_gamma(i) = (alpha * diffA(t_u)) * t_grad(i);
625 // take fist element to local matrix
627 &*locMat.data().begin());
628 // loop over rows base functions
629 for (int rr = 0; rr != nbRows; rr++) {
630 // get column base function
631 auto t_col = col_data.getFTensor0N(gg, 0);
632 // get column base functions gradient at gauss point gg
633 auto t_col_grad = col_data.getFTensor1DiffN<3>(gg, 0);
634 // loop over columns
635 for (int cc = 0; cc != nbCols; cc++) {
636 // calculate element of local matrix
637 a += (t_row_grad(i) * beta) * t_col_grad(i) +
638 t_row_grad(i) * (t_gamma(i) * t_col);
639 ++t_col; // move to next base function
640 ++t_col_grad; // move to another gradient of base function on column
641 ++a; // move to another element of local matrix in column
642 }
643 ++t_row_grad; // move to another element of gradient of base function on
644 // row
645 }
646 ++t_w; // move to another integration weight
647 ++t_u; // move to next value at integration point
648 ++t_grad; // move to next gradient value
649 }
651 }
652
653 boost::function<double(const double)> A;
654 boost::function<double(const double)> diffA;
655 boost::shared_ptr<VectorDouble> fieldVals;
656 boost::shared_ptr<MatrixDouble> gradVals;
657};
658
659struct OpResF_Domain : public OpF {
660
661 OpResF_Domain(FSource f_source, boost::function<double(const double)> a,
662 boost::shared_ptr<VectorDouble> &field_vals,
663 boost::shared_ptr<MatrixDouble> &grad_vals)
664 : OpF(f_source), A(a), fieldVals(field_vals), gradVals(grad_vals) {}
665
666protected:
668
669 /**
670 * \brief Integrate local entity vector
671 * @param data entity data on element row
672 * @return error code
673 */
674 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
676 // set size of local vector
677 locVec.resize(nbRows, false);
678 // clear local entity vector
679 locVec.clear();
680 // get finite element volume
681 double vol = getVolume();
682 // get integration weights
683 auto t_w = getFTensor0IntegrationWeight();
684 // get solution at integration point
685 auto t_u = getFTensor0FromVec(*fieldVals);
686 // get solution at integration point
687 auto t_grad = getFTensor1FromMat<3>(*gradVals);
688 // get base functions on entity
689 auto t_v = data.getFTensor0N();
690 // get base function gradient on rows
691 auto t_v_grad = data.getFTensor1DiffN<3>();
692 // get coordinates at integration points
693 auto t_coords = getFTensor1CoordsAtGaussPts();
694 // loop over all integration points
695 for (int gg = 0; gg != nbIntegrationPts; gg++) {
696 // evaluate constant term
697 const double alpha = vol * t_w;
698 const double source_term =
699 alpha * fSource(t_coords(NX), t_coords(NY), t_coords(NZ));
701 grad_term(i) = (alpha * A(t_u)) * t_grad(i);
702 // get element of local vector
704 &*locVec.data().begin());
705 // loop over base functions
706 for (int rr = 0; rr != nbRows; rr++) {
707 // add to local vector source term
708 t_a += t_v_grad(i) * grad_term(i) + t_v * source_term;
709 ++t_a; // move to next element of local vector
710 ++t_v; // move to next base function
711 ++t_v_grad; // move to next gradient of base function
712 }
713 ++t_w; // move to next integration weights
714 ++t_u; // move to next value
715 ++t_grad; // move to next gradient value
716 ++t_coords; // move to next physical coordinates at integration point
717 }
719 }
720
721 boost::function<double(const double)> A;
722 boost::shared_ptr<VectorDouble> fieldVals;
723 boost::shared_ptr<MatrixDouble> gradVals;
724};
725
726struct OpRes_g : public Op_g {
727
728 OpRes_g(FVal f_value, boost::shared_ptr<VectorDouble> &field_vals)
729 : Op_g(f_value, "L", 1), fieldVals(field_vals) {}
730
731protected:
732 boost::shared_ptr<VectorDouble> fieldVals;
733
734 /**
735 * \brief Integrate local constrains vector
736 */
737 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
739 // set size to local vector
740 locVec.resize(nbRows, false);
741 // clear local vector
742 locVec.clear();
743 // get face area
744 const double area = getArea() * bEta;
745 // get integration weight
746 auto t_w = getFTensor0IntegrationWeight();
747 // get base function
748 auto t_l = data.getFTensor0N();
749 // get solution at integration point
750 auto t_u = getFTensor0FromVec(*fieldVals);
751 // get coordinates at integration point
752 auto t_coords = getFTensor1CoordsAtGaussPts();
753 // make loop over integration points
754 for (int gg = 0; gg != nbIntegrationPts; gg++) {
755 // evaluate function on boundary and scale it by area and integration
756 // weight
757 double alpha = area * t_w;
758 // get element of vector
760 &*locVec.data().begin());
761 for (int rr = 0; rr != nbRows; rr++) {
762 t_a += alpha * t_l *
763 (t_u - fValue(t_coords(NX), t_coords(NY), t_coords(NZ)));
764 ++t_a;
765 ++t_l;
766 }
767 ++t_w;
768 ++t_u;
769 ++t_coords;
770 }
772 }
773};
774
775struct OpResF_Boundary : public Op_g {
776
777 OpResF_Boundary(boost::shared_ptr<VectorDouble> &lambda_vals)
778 : Op_g(FVal(), "U", 1), lambdaVals(lambda_vals) {}
779
780protected:
781 boost::shared_ptr<VectorDouble> lambdaVals;
782
783 /**
784 * \brief Integrate local constrains vector
785 */
786 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
788 // set size to local vector
789 locVec.resize(nbRows, false);
790 // clear local vector
791 locVec.clear();
792 // get face area
793 const double area = getArea() * bEta;
794 // get integration weight
795 auto t_w = getFTensor0IntegrationWeight();
796 // get base function
797 auto t_u = data.getFTensor0N();
798 // get solution at integration point
799 auto t_lambda = getFTensor0FromVec(*lambdaVals);
800 // make loop over integration points
801 for (int gg = 0; gg != nbIntegrationPts; gg++) {
802 // evaluate function on boundary and scale it by area and integration
803 // weight
804 double alpha = area * t_w;
805 // get element of vector
807 &*locVec.data().begin());
808 for (int rr = 0; rr != nbRows; rr++) {
809 t_a += alpha * t_u * t_lambda;
810 ++t_a;
811 ++t_u;
812 }
813 ++t_w;
814 ++t_lambda;
815 }
817 }
818};
819
820/**
821 * \brief Set integration rule to volume elements
822 *
823 * This rule is used to integrate \f$\nabla v \cdot \nabla u\f$, thus
824 * if approximation field and testing field are polynomial order "p",
825 * then rule for exact integration is 2*(p-1).
826 *
827 * Integration rule is order of polynomial which is calculated exactly. Finite
828 * element selects integration method based on return of this function.
829 *
830 */
831struct VolRule {
832 int operator()(int, int, int p) const { return 2 * (p - 1); }
833};
834
835/**
836 * \brief Set integration rule to boundary elements
837 *
838 * This is uses to integrate values on the face. Is used to integrate
839 * \f$(\mathbf{n} \cdot \lambda) u\f$, where Lagrange multiplayer
840 * is order "p_row" and approximate function is order "p_col".
841 *
842 * Integration rule is order of polynomial which is calculated exactly. Finite
843 * element selects integration method based on return of this function.
844 *
845 */
846struct FaceRule {
847 int operator()(int p_row, int p_col, int p_data) const {
848 return 2 * p_data + 1;
849 }
850};
851
852/**
853 * \brief Create finite elements instances
854 *
855 * Create finite element instances and add operators to finite elements.
856 *
857 */
859
861
862 /**
863 * \brief Create finite element to calculate matrix and vectors
864 */
866 boost::function<double(const double, const double, const double)> f_u,
867 boost::function<double(const double, const double, const double)>
868 f_source,
869 boost::shared_ptr<ForcesAndSourcesCore> &domain_lhs_fe,
870 boost::shared_ptr<ForcesAndSourcesCore> &boundary_lhs_fe,
871 boost::shared_ptr<ForcesAndSourcesCore> &domain_rhs_fe,
872 boost::shared_ptr<ForcesAndSourcesCore> &boundary_rhs_fe,
873 bool trans = true) const {
875
876 // Create elements element instances
877 domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
878 new VolumeElementForcesAndSourcesCore(mField));
879 boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
881 domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
882 new VolumeElementForcesAndSourcesCore(mField));
883 boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
885
886 // Set integration rule to elements instances
887 domain_lhs_fe->getRuleHook = VolRule();
888 domain_rhs_fe->getRuleHook = VolRule();
889 boundary_lhs_fe->getRuleHook = FaceRule();
890 boundary_rhs_fe->getRuleHook = FaceRule();
891
892 // Add operators to element instances
893 // Add operator grad-grad for calculate matrix
894 domain_lhs_fe->getOpPtrVector().push_back(new OpK());
895 // Add operator to calculate source terms
896 domain_rhs_fe->getOpPtrVector().push_back(new OpF(f_source));
897 // Add operator calculating constrains matrix
898 boundary_lhs_fe->getOpPtrVector().push_back(new OpC(trans));
899 // Add operator calculating constrains vector
900 boundary_rhs_fe->getOpPtrVector().push_back(new Op_g(f_u));
901
903 }
904
905 /**
906 * \brief Create finite element to calculate error
907 */
909 boost::function<double(const double, const double, const double)> f_u,
910 boost::function<FTensor::Tensor1<double, 3>(const double, const double,
911 const double)>
912 g_u,
913 Vec global_error,
914 boost::shared_ptr<ForcesAndSourcesCore> &domain_error) const {
916 // Create finite element instance to calculate error
917 domain_error = boost::shared_ptr<ForcesAndSourcesCore>(
918 new VolumeElementForcesAndSourcesCore(mField));
919 domain_error->getRuleHook = VolRule();
920 // Set integration rule
921 // Crate shared vector storing values of field "u" on integration points on
922 // element. element is local and is used to exchange data between operators.
923 boost::shared_ptr<VectorDouble> values_at_integration_ptr =
924 boost::make_shared<VectorDouble>();
925 // Storing gradients of field
926 boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
927 boost::make_shared<MatrixDouble>();
928 // Add default operator to calculate field values at integration points
929 domain_error->getOpPtrVector().push_back(
930 new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
931 // Add default operator to calculate field gradient at integration points
932 domain_error->getOpPtrVector().push_back(
933 new OpCalculateScalarFieldGradient<3>("U", grad_at_integration_ptr));
934 // Add operator to integrate error element by element.
935 domain_error->getOpPtrVector().push_back(
936 new OpError(f_u, g_u, values_at_integration_ptr,
937 grad_at_integration_ptr, global_error));
939 }
940
941 /**
942 * \brief Create finite element to post-process results
943 */
945 boost::shared_ptr<PostProcFE> &post_proc_volume) const {
946
948
949 // Note that user can stack together arbitrary number of operators to
950 // compose complex PDEs.
951
952 // Post-process results. This is standard element, with functionality
953 // enabling refining mesh for post-processing. In addition in
954 // PostProcOnRefMesh.hpp are implanted set of users operators to
955 // post-processing fields. Here using simplified mechanism for
956 // post-processing finite element, we add operators to save data from field
957 // on mesh tags for ParaView visualization.
958 post_proc_volume =
959 boost::shared_ptr<PostProcFE>(new PostProcFE(mField));
960 // Add operators to the elements, starting with some generic
961
962 auto det_ptr = boost::make_shared<VectorDouble>();
963 auto jac_ptr = boost::make_shared<MatrixDouble>();
964 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
965 post_proc_volume->getOpPtrVector().push_back(
966 new OpCalculateHOJac<3>(jac_ptr));
967 post_proc_volume->getOpPtrVector().push_back(
968 new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
969 post_proc_volume->getOpPtrVector().push_back(
970 new OpSetHOInvJacToScalarBases<3>(H1, inv_jac_ptr));
971
972 auto u_ptr = boost::make_shared<VectorDouble>();
973 auto grad_ptr = boost::make_shared<MatrixDouble>();
974 auto e_ptr = boost::make_shared<VectorDouble>();
975
976 post_proc_volume->getOpPtrVector().push_back(
977 new OpCalculateScalarFieldValues("U", u_ptr));
978 post_proc_volume->getOpPtrVector().push_back(
979 new OpCalculateScalarFieldGradient<3>("U", grad_ptr));
980 post_proc_volume->getOpPtrVector().push_back(
981 new OpCalculateScalarFieldValues("ERROR", e_ptr));
982
983 using OpPPMap = OpPostProcMapInMoab<3, 3>;
984
985 post_proc_volume->getOpPtrVector().push_back(
986
987 new OpPPMap(
988
989 post_proc_volume->getPostProcMesh(),
990 post_proc_volume->getMapGaussPts(),
991
992 {{"U", u_ptr}, {"ERROR", e_ptr}},
993
994 {{"GRAD", grad_ptr}},
995
996 {},
997
998 {}
999
1000 )
1001
1002 );
1003
1005 }
1006
1007 /**
1008 * \brief Create finite element to calculate matrix and vectors
1009 */
1011 boost::function<double(const double, const double, const double)> f_u,
1012 boost::function<double(const double, const double, const double)>
1013 f_source,
1014 boost::function<double(const double)> a,
1015 boost::function<double(const double)> diff_a,
1016 boost::shared_ptr<ForcesAndSourcesCore> &domain_lhs_fe,
1017 boost::shared_ptr<ForcesAndSourcesCore> &boundary_lhs_fe,
1018 boost::shared_ptr<ForcesAndSourcesCore> &domain_rhs_fe,
1019 boost::shared_ptr<ForcesAndSourcesCore> &boundary_rhs_fe,
1020 ForcesAndSourcesCore::RuleHookFun vol_rule,
1021 ForcesAndSourcesCore::RuleHookFun face_rule = FaceRule(),
1022 bool trans = true) const {
1024
1025 // Create elements element instances
1026 domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1027 new VolumeElementForcesAndSourcesCore(mField));
1028 boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1030 domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1031 new VolumeElementForcesAndSourcesCore(mField));
1032 boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1034
1035 // Set integration rule to elements instances
1036 domain_lhs_fe->getRuleHook = vol_rule;
1037 domain_rhs_fe->getRuleHook = vol_rule;
1038 boundary_lhs_fe->getRuleHook = face_rule;
1039 boundary_rhs_fe->getRuleHook = face_rule;
1040
1041 // Set integration rule
1042 // Crate shared vector storing values of field "u" on integration points on
1043 // element. element is local and is used to exchange data between operators.
1044 boost::shared_ptr<VectorDouble> values_at_integration_ptr =
1045 boost::make_shared<VectorDouble>();
1046 // Storing gradients of field
1047 boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
1048 boost::make_shared<MatrixDouble>();
1049 // multipliers values
1050 boost::shared_ptr<VectorDouble> multiplier_at_integration_ptr =
1051 boost::make_shared<VectorDouble>();
1052
1053 // Add operators to element instances
1054 // Add default operator to calculate field values at integration points
1055 domain_lhs_fe->getOpPtrVector().push_back(
1056 new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
1057 // Add default operator to calculate field gradient at integration points
1058 domain_lhs_fe->getOpPtrVector().push_back(
1059 new OpCalculateScalarFieldGradient<3>("U", grad_at_integration_ptr));
1060 // Add operator grad-(1+u^2)grad for calculate matrix
1061 domain_lhs_fe->getOpPtrVector().push_back(new OpKt(
1062 a, diff_a, values_at_integration_ptr, grad_at_integration_ptr));
1063
1064 // Add default operator to calculate field values at integration points
1065 domain_rhs_fe->getOpPtrVector().push_back(
1066 new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
1067 // Add default operator to calculate field gradient at integration points
1068 domain_rhs_fe->getOpPtrVector().push_back(
1069 new OpCalculateScalarFieldGradient<3>("U", grad_at_integration_ptr));
1070 // Add operator to calculate source terms
1071 domain_rhs_fe->getOpPtrVector().push_back(new OpResF_Domain(
1072 f_source, a, values_at_integration_ptr, grad_at_integration_ptr));
1073
1074 // Add operator calculating constrains matrix
1075 boundary_lhs_fe->getOpPtrVector().push_back(new OpC(trans));
1076
1077 // Add default operator to calculate field values at integration points
1078 boundary_rhs_fe->getOpPtrVector().push_back(
1079 new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
1080 // Add default operator to calculate values of Lagrange multipliers
1081 boundary_rhs_fe->getOpPtrVector().push_back(
1082 new OpCalculateScalarFieldValues("L", multiplier_at_integration_ptr));
1083 // Add operator calculating constrains vector
1084 boundary_rhs_fe->getOpPtrVector().push_back(
1085 new OpRes_g(f_u, values_at_integration_ptr));
1086 boundary_rhs_fe->getOpPtrVector().push_back(
1087 new OpResF_Boundary(multiplier_at_integration_ptr));
1088
1090 }
1091
1092private:
1094};
1095
1096} // namespace PoissonExample
1097
1098#endif //__POISSONOPERATORS_HPP__
static Index< 'p', 3 > p
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ 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
const double c
speed of light (cm/ns)
PostProcBrokenMeshInMoab< VolumeElementForcesAndSourcesCore > PostProcFE
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto field_name
Deprecated interface functions.
Post post-proc data at points from hash maps.
Create finite elements instances.
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< PostProcFE > &post_proc_volume) const
Create finite element to post-process results.
CreateFiniteElements(MoFEM::Interface &m_field)
MoFEMErrorCode createFEToAssembleMatrixAndVectorForNonlinearProblem(boost::function< double(const double, const double, const double)> f_u, boost::function< double(const double, const double, const double)> f_source, boost::function< double(const double)> a, boost::function< double(const double)> diff_a, boost::shared_ptr< ForcesAndSourcesCore > &domain_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &domain_rhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_rhs_fe, ForcesAndSourcesCore::RuleHookFun vol_rule, ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule(), bool trans=true) const
Create finite element to calculate matrix and vectors.
MoFEMErrorCode createFEToEvaluateError(boost::function< double(const double, const double, const double)> f_u, boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> g_u, Vec global_error, boost::shared_ptr< ForcesAndSourcesCore > &domain_error) const
Create finite element to calculate error.
MoFEMErrorCode createFEToAssembleMatrixAndVector(boost::function< double(const double, const double, const double)> f_u, boost::function< double(const double, const double, const double)> f_source, boost::shared_ptr< ForcesAndSourcesCore > &domain_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &domain_rhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_rhs_fe, bool trans=true) const
Create finite element to calculate matrix and vectors.
Set integration rule to boundary elements.
int operator()(int p_row, int p_col, int p_data) const
Assemble constrains vector.
FTensor::Number< 2 > NZ
z-direction index
Op_g(FVal f_value, const string field_name="L", const double beta=1)
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
assemble constrains vectors
boost::function< double(const double, const double, const double)> FVal
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
FTensor::Number< 0 > NX
x-direction index
FTensor::Number< 1 > NY
y-direction index
FVal fValue
Function pointer evaluating values of "U" at the boundary.
template class for integration oh the right hand side
int nbIntegrationPts
number of integration points
virtual MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)=0
Class dedicated to assemble operator to global system vector.
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
This function is called by finite element.
OpBaseRhs(const std::string field_name)
virtual MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)=0
Class dedicated to integrate operator.
Calculate constrains matrix.
MatrixDouble locMat
local constrains matrix
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
integrate local constrains matrix
OpC(const bool assemble_transpose)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate local constrains matrix.
int nbIntegrationPts
number of integration points
int nbRows
< error code
const bool assembleTranspose
assemble transpose, i.e. CT if set to true
int nbCols
number of dofs on column
Vec globalError
ghost vector with global (integrated over volume) error
UVal uValue
function with exact solution
OpError(UVal u_value, GVal g_value, boost::shared_ptr< VectorDouble > &field_vals, boost::shared_ptr< MatrixDouble > &grad_vals, Vec global_error)
boost::shared_ptr< MatrixDouble > gradVals
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
boost::function< double(const double, const double, const double)> UVal
FTensor::Number< 2 > NZ
FTensor::Index< 'i', 3 > i
GVal gValue
function with exact solution for gradient
boost::shared_ptr< VectorDouble > fieldVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate error.
FTensor::Number< 1 > NY
boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> GVal
FTensor::Number< 0 > NX
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
Assemble error.
Operator calculate source term,.
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
assemble local entity vector to the global right hand side
FTensor::Number< 2 > NZ
OpF(FSource f_source)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local entity vector.
boost::function< double(const double, const double, const double)> FSource
FTensor::Number< 1 > NY
FTensor::Number< 0 > NX
Calculate the grad-grad operator and assemble matrix.
virtual MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
int nbRows
< error code
int nbCols
number if dof on column
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations for give operator.
bool isDiag
true if this block is on diagonal
MatrixDouble locMat
local entity block matrix
int nbIntegrationPts
number of integration points
FTensor::Index< 'i', 3 > i
summit Index
virtual MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble local entity block matrix.
boost::function< double(const double)> diffA
boost::shared_ptr< MatrixDouble > gradVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< VectorDouble > fieldVals
boost::function< double(const double)> A
OpKt(boost::function< double(const double)> a, boost::function< double(const double)> diff_a, boost::shared_ptr< VectorDouble > &field_vals, boost::shared_ptr< MatrixDouble > &grad_vals)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
boost::shared_ptr< VectorDouble > fieldVals
OpRes_g(FVal f_value, boost::shared_ptr< VectorDouble > &field_vals)
boost::shared_ptr< VectorDouble > lambdaVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
OpResF_Boundary(boost::shared_ptr< VectorDouble > &lambda_vals)
boost::shared_ptr< MatrixDouble > gradVals
boost::shared_ptr< VectorDouble > fieldVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local entity vector.
boost::function< double(const double)> A
FTensor::Index< 'i', 3 > i
OpResF_Domain(FSource f_source, boost::function< double(const double)> a, boost::shared_ptr< VectorDouble > &field_vals, boost::shared_ptr< MatrixDouble > &grad_vals)
Set integration rule to volume elements.
int operator()(int, int, int p) const