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