v0.14.0
Loading...
Searching...
No Matches
Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | List of all members
OpK Struct Reference
Inheritance diagram for OpK:
[legend]
Collaboration diagram for OpK:
[legend]

Public Member Functions

 OpK (bool symm=true)
 
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. More...
 

Public Attributes

MatrixDouble K
 
FTensor::Ddg< double, 3, 3 > tD
 
double yOung
 
double pOisson
 

Protected Member Functions

MoFEMErrorCode iNtegrate (EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
 Integrate B^T D B operator. More...
 
MoFEMErrorCode aSsemble (EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
 Assemble local entity block matrix. More...
 

Protected Attributes

int nbRows
 number of dofs on rows More...
 
int nbCols
 number if dof on column More...
 
int nbIntegrationPts
 number of integration points More...
 
bool isDiag
 true if this block is on diagonal More...
 

Detailed Description

Examples
PoissonOperators.hpp, eigen_elastic.cpp, and simple_elasticity.cpp.

Definition at line 16 of file simple_elasticity.cpp.

Constructor & Destructor Documentation

◆ OpK()

OpK::OpK ( bool  symm = true)
inline

Definition at line 29 of file simple_elasticity.cpp.

30 : VolumeElementForcesAndSourcesCore::UserDataOperator("U", "U", OPROWCOL,
31 symm) {
32
33 // Evaluation of the elastic stiffness tensor, D
34
35 // hardcoded choice of elastic parameters
36 pOisson = 0.1;
37 yOung = 10;
38
39 // coefficient used in intermediate calculation
40 const double coefficient = yOung / ((1 + pOisson) * (1 - 2 * pOisson));
41
46
47 tD(i, j, k, l) = 0.;
48
49 tD(0, 0, 0, 0) = 1 - pOisson;
50 tD(1, 1, 1, 1) = 1 - pOisson;
51 tD(2, 2, 2, 2) = 1 - pOisson;
52
53 tD(0, 1, 0, 1) = 0.5 * (1 - 2 * pOisson);
54 tD(0, 2, 0, 2) = 0.5 * (1 - 2 * pOisson);
55 tD(1, 2, 1, 2) = 0.5 * (1 - 2 * pOisson);
56
57 tD(0, 0, 1, 1) = pOisson;
58 tD(1, 1, 0, 0) = pOisson;
59 tD(0, 0, 2, 2) = pOisson;
60 tD(2, 2, 0, 0) = pOisson;
61 tD(1, 1, 2, 2) = pOisson;
62 tD(2, 2, 1, 1) = pOisson;
63
64 tD(i, j, k, l) *= coefficient;
65 }
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double yOung
double pOisson
FTensor::Ddg< double, 3, 3 > tD

Member Function Documentation

◆ aSsemble()

MoFEMErrorCode OpK::aSsemble ( EntitiesFieldData::EntData row_data,
EntitiesFieldData::EntData col_data 
)
inlineprotected

Assemble local entity block matrix.

Parameters
row_datarow data (consist base functions on row entity)
col_datacolumn data (consist base functions on column entity)
Returns
error code
Examples
simple_elasticity.cpp.

Definition at line 205 of file simple_elasticity.cpp.

206 {
208 // get pointer to first global index on row
209 const int *row_indices = &*row_data.getIndices().data().begin();
210 // get pointer to first global index on column
211 const int *col_indices = &*col_data.getIndices().data().begin();
212 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
213 : getFEMethod()->snes_B;
214 // assemble local matrix
215 CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
216 &*K.data().begin(), ADD_VALUES);
217
218 if (!isDiag && sYmm) {
219 // if not diagonal term and since global matrix is symmetric assemble
220 // transpose term.
221 K = trans(K);
222 CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
223 &*K.data().begin(), ADD_VALUES);
224 }
226 }
#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
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
bool isDiag
true if this block is on diagonal
int nbRows
number of dofs on rows
MatrixDouble K
int nbCols
number if dof on column

◆ doWork()

MoFEMErrorCode OpK::doWork ( int  row_side,
int  col_side,
EntityType  row_type,
EntityType  col_type,
EntitiesFieldData::EntData row_data,
EntitiesFieldData::EntData col_data 
)
inline

Do calculations for give operator.

Parameters
row_siderow side number (local number) of entity on element
col_sidecolumn side number (local number) of entity on element
row_typetype of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
col_typetype of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
row_datadata for row
col_datadata for column
Returns
error code
Examples
simple_elasticity.cpp.

Definition at line 77 of file simple_elasticity.cpp.

80 {
81
83
84 // get number of dofs on row
85 nbRows = row_data.getIndices().size();
86 // if no dofs on row, exit that work, nothing to do here
87 if (!nbRows)
89
90 // get number of dofs on column
91 nbCols = col_data.getIndices().size();
92 // if no dofs on Columbia, exit nothing to do here
93 if (!nbCols)
95
96 // K_ij matrix will have 3 times the number of degrees of freedom of the
97 // i-th entity set (nbRows)
98 // and 3 times the number of degrees of freedom of the j-th entity set
99 // (nbCols)
100 K.resize(nbRows, nbCols, false);
101 K.clear();
102
103 // get number of integration points
104 nbIntegrationPts = getGaussPts().size2();
105 // check if entity block is on matrix diagonal
106 if (row_side == col_side && row_type == col_type) {
107 isDiag = true;
108 } else {
109 isDiag = false;
110 }
111
112 // integrate local matrix for entity block
113 CHKERR iNtegrate(row_data, col_data);
114
115 // assemble local matrix
116 CHKERR aSsemble(row_data, col_data);
117
119 }
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate B^T D B operator.
int nbIntegrationPts
number of integration points
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble local entity block matrix.

◆ iNtegrate()

MoFEMErrorCode OpK::iNtegrate ( EntitiesFieldData::EntData row_data,
EntitiesFieldData::EntData col_data 
)
inlineprotected

Integrate B^T D B operator.

Parameters
row_datarow data (consist base functions on row entity)
col_datacolumn data (consist base functions on column entity)
Returns
error code
Examples
simple_elasticity.cpp.

Definition at line 133 of file simple_elasticity.cpp.

134 {
136
137 // get sub-block (3x3) of local stiffens matrix, here represented by second
138 // order tensor
139 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
141 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
142 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
143 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
144 };
145
150
151 // get element volume
152 double vol = getVolume();
153
154 // get intergration weights
155 auto t_w = getFTensor0IntegrationWeight();
156
157 // get derivatives of base functions on rows
158 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
159 // iterate over integration points
160 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
161
162 // calculate scalar weight times element volume
163 const double a = t_w * vol;
164
165 // iterate over row base functions
166 for (int rr = 0; rr != nbRows / 3; ++rr) {
167
168 // get sub matrix for the row
169 auto t_m = get_tensor2(K, 3 * rr, 0);
170
171 // get derivatives of base functions for columns
172 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
173
174 // iterate column base functions
175 for (int cc = 0; cc != nbCols / 3; ++cc) {
176
177 // integrate block local stiffens matrix
178 t_m(i, k) +=
179 a * (tD(i, j, k, l) * (t_row_diff_base(j) * t_col_diff_base(l)));
180
181 // move to next column base function
182 ++t_col_diff_base;
183
184 // move to next block of local stiffens matrix
185 ++t_m;
186 }
187
188 // move to next row base function
189 ++t_row_diff_base;
190 }
191
192 // move to next integration weight
193 ++t_w;
194 }
195
197 }
constexpr double a
FTensor::Index< 'm', SPACE_DIM > m
const double c
speed of light (cm/ns)
int r
Definition: sdf.py:8
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.

Member Data Documentation

◆ isDiag

bool OpK::isDiag
protected

true if this block is on diagonal

Examples
simple_elasticity.cpp.

Definition at line 125 of file simple_elasticity.cpp.

◆ K

MatrixDouble OpK::K
Examples
simple_elasticity.cpp.

Definition at line 19 of file simple_elasticity.cpp.

◆ nbCols

int OpK::nbCols
protected

number if dof on column

Examples
simple_elasticity.cpp.

Definition at line 123 of file simple_elasticity.cpp.

◆ nbIntegrationPts

int OpK::nbIntegrationPts
protected

number of integration points

Examples
simple_elasticity.cpp.

Definition at line 124 of file simple_elasticity.cpp.

◆ nbRows

int OpK::nbRows
protected

number of dofs on rows

Examples
simple_elasticity.cpp.

Definition at line 122 of file simple_elasticity.cpp.

◆ pOisson

double OpK::pOisson
Examples
simple_elasticity.cpp.

Definition at line 27 of file simple_elasticity.cpp.

◆ tD

FTensor::Ddg<double, 3, 3> OpK::tD
Examples
simple_elasticity.cpp.

Definition at line 22 of file simple_elasticity.cpp.

◆ yOung

double OpK::yOung
Examples
simple_elasticity.cpp.

Definition at line 25 of file simple_elasticity.cpp.


The documentation for this struct was generated from the following file: