v0.13.1
Loading...
Searching...
No Matches
PoissonOperators.hpp
Go to the documentation of this file.
1#ifndef __POISSONOPERTATORS_HPP__
2#define __POISSONOPERTATORS_HPP__
3
4#include <stdlib.h>
6
7namespace PoissonOps {
8
10using OpVolEle = FaceElementForcesAndSourcesCore::UserDataOperator;
11
14
15using EntData = DataForcesAndSourcesCore::EntData;
16
17// Operators
18
20
21const double body_source = 1;
22
23struct OpLhsUU : public OpVolEle {
24public:
25 OpLhsUU(std::string domain_field1, std::string domain_field2)
26 : OpVolEle(domain_field1, domain_field2, OpVolEle::OPROWCOL) {
27 sYmm = true;
28 }
29
30 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
31 EntityType col_type, EntData &row_data,
32 EntData &col_data) {
34
35 const int nb_row_dofs = row_data.getIndices().size();
36 const int nb_col_dofs = col_data.getIndices().size();
37
38 if (nb_row_dofs && nb_col_dofs) {
39 locLhs.resize(nb_row_dofs, nb_col_dofs, false);
40 locLhs.clear();
41 const int nb_integration_pts = getGaussPts().size2();
43
44 auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
45
46 const double vol = getMeasure();
47 for (int gg = 0; gg != nb_integration_pts; ++gg) {
48 const double a = t_w * vol;
49 for (int rr = 0; rr != nb_row_dofs; ++rr) {
50 auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
51 for (int cc = 0; cc != nb_col_dofs; ++cc) {
52 locLhs(rr, cc) += a * t_row_diff_base(i) * t_col_diff_base(i);
53 ++t_col_diff_base;
54 } // endFor cc
55 ++t_row_diff_base;
56 } // endFor rr
57 ++t_w;
58 } // endFor gg
59 CHKERR MatSetValues(getFEMethod()->ksp_B, row_data, col_data,
60 &locLhs(0, 0), ADD_VALUES);
61 if (row_side != col_side || row_type != col_type) {
62 transLocLhs.resize(nb_col_dofs, nb_row_dofs, false);
63 noalias(transLocLhs) = trans(locLhs);
64 CHKERR MatSetValues(getFEMethod()->ksp_B, col_data, row_data,
65 &transLocLhs(0, 0), ADD_VALUES);
66 } // (row_side != col_side || row_type != col_type)
67 } // endIf (nb_row_dofs && nb_col_dofs)
68
70 }
71
72private:
74};
75
76struct OpRhsU : public OpVolEle {
77public:
78 OpRhsU(std::string domain_field) : OpVolEle(domain_field, OpVolEle::OPROW) {}
79 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
81
82 const int nb_dofs = data.getIndices().size();
83 if (nb_dofs) {
84 locRhs.resize(nb_dofs, false);
85 locRhs.clear();
86 const int nb_integration_pts = getGaussPts().size2();
87
88 auto t_base = data.getFTensor0N();
90 const double vol = getMeasure();
91 for (int gg = 0; gg != nb_integration_pts; ++gg) {
92 const double a = vol * t_w;
93
94 for (int rr = 0; rr != nb_dofs; ++rr) {
95 locRhs[rr] += a * (t_base * body_source);
96 ++t_base;
97 } // endFor rr
98 ++t_w;
99 } // endFor gg
100 CHKERR VecSetOption(getFEMethod()->ksp_f, VEC_IGNORE_NEGATIVE_INDICES,
101 PETSC_TRUE);
102 CHKERR VecSetValues(getFEMethod()->ksp_f, data, &*locRhs.begin(),
103 ADD_VALUES);
104 } // endIf (nb_dofs)
105
107 }
108
109private:
111};
112
113}; // namespace PoissonOps
114
115#endif //__POISSONOPERTATORS_HPP__
constexpr double a
#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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
DataForcesAndSourcesCore::EntData EntData
FTensor::Index< 'i', 2 > i
const double body_source
bool sYmm
If true assume that matrix is symmetric structure.
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpLhsUU(std::string domain_field1, std::string domain_field2)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpRhsU(std::string domain_field)