v0.9.0
Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | List of all members
OpK Struct Reference

assemble the stiffness matrix More...

#include <users_modules/linear_isotropic_hardening/src/LinearIsotropicHardening.hpp>

Inherits UserDataOperator, and UserDataOperator.

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, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
 Do calculations for give operator. More...
 
 OpK (CommonData &common_data, bool symm=true)
 
MoFEMErrorCode doWork (int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
 Do calculations for give operator. More...
 

Public Attributes

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

Protected Member Functions

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

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

assemble the stiffness matrix

Examples
PoissonOperators.hpp, and simple_elasticity.cpp.

Definition at line 29 of file simple_elasticity.cpp.

Constructor & Destructor Documentation

◆ OpK() [1/2]

OpK::OpK ( bool  symm = true)

Definition at line 42 of file simple_elasticity.cpp.

44  symm) {
45 
46  // Evaluation of the elastic stiffness tensor, D
47 
48  // hardcoded choice of elastic parameters
49  pOisson = 0.1;
50  yOung = 10;
51 
52  // coefficient used in intermediate calculation
53  const double coefficient = yOung / ((1 + pOisson) * (1 - 2 * pOisson));
54 
59 
60  tD(i, j, k, l) = 0.;
61 
62  tD(0, 0, 0, 0) = 1 - pOisson;
63  tD(1, 1, 1, 1) = 1 - pOisson;
64  tD(2, 2, 2, 2) = 1 - pOisson;
65 
66  tD(0, 1, 0, 1) = 0.5 * (1 - 2 * pOisson);
67  tD(0, 2, 0, 2) = 0.5 * (1 - 2 * pOisson);
68  tD(1, 2, 1, 2) = 0.5 * (1 - 2 * pOisson);
69 
70  tD(0, 0, 1, 1) = pOisson;
71  tD(1, 1, 0, 0) = pOisson;
72  tD(0, 0, 2, 2) = pOisson;
73  tD(2, 2, 0, 0) = pOisson;
74  tD(1, 1, 2, 2) = pOisson;
75  tD(2, 2, 1, 1) = pOisson;
76 
77  tD(i, j, k, l) *= coefficient;
78 
79  }
double pOisson
FTensor::Ddg< double, 3, 3 > tD
double yOung
ForcesAndSourcesCore::UserDataOperator UserDataOperator

◆ OpK() [2/2]

OpK::OpK ( CommonData common_data,
bool  symm = true 
)

Definition at line 66 of file LinearIsotropicHardening.hpp.

68  symm),
69  commonData(common_data) {}
CommonData & commonData
ForcesAndSourcesCore::UserDataOperator UserDataOperator

Member Function Documentation

◆ aSsemble() [1/2]

MoFEMErrorCode OpK::aSsemble ( DataForcesAndSourcesCore::EntData row_data,
DataForcesAndSourcesCore::EntData col_data 
)
protected

Definition at line 202 of file LinearIsotropicHardening.hpp.

203  {
205  // get pointer to first global index on row
206  const int *row_indices = &*row_data.getIndices().data().begin();
207  // get pointer to first global index on column
208  const int *col_indices = &*col_data.getIndices().data().begin();
209  // Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
210  // : getFEMethod()->snes_B;
211  // assemble local matrix
212  CHKERR MatSetValues(getFEMethod()->snes_B, nbRows, row_indices, nbCols, col_indices,
213  &*K.data().begin(), ADD_VALUES);
214 
215  if (!isDiag && sYmm) {
216  // if not diagonal term and since global matrix is symmetric assemble
217  // transpose term.
218  K = trans(K);
219  CHKERR MatSetValues(getFEMethod()->snes_B, nbCols, col_indices, nbRows, row_indices,
220  &*K.data().begin(), ADD_VALUES);
221  }
223  }
bool isDiag
true if this block is on diagonal
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MatrixDouble K
int nbRows
number of dofs on rows
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
#define CHKERR
Inline error check.
Definition: definitions.h:596
int nbCols
number if dof on column
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ aSsemble() [2/2]

MoFEMErrorCode OpK::aSsemble ( DataForcesAndSourcesCore::EntData row_data,
DataForcesAndSourcesCore::EntData col_data 
)
protected

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

Definition at line 220 of file simple_elasticity.cpp.

221  {
223  // get pointer to first global index on row
224  const int *row_indices = &*row_data.getIndices().data().begin();
225  // get pointer to first global index on column
226  const int *col_indices = &*col_data.getIndices().data().begin();
227  Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
228  : getFEMethod()->snes_B;
229  // assemble local matrix
230  CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
231  &*K.data().begin(), ADD_VALUES);
232 
233  if (!isDiag && sYmm) {
234  // if not diagonal term and since global matrix is symmetric assemble
235  // transpose term.
236  K = trans(K);
237  CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
238  &*K.data().begin(), ADD_VALUES);
239  }
241  }
bool isDiag
true if this block is on diagonal
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MatrixDouble K
const VectorInt & getIndices() const
Get global indices of dofs on entity.
int nbRows
number of dofs on rows
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
#define CHKERR
Inline error check.
Definition: definitions.h:596
int nbCols
number if dof on column
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ doWork() [1/2]

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

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

Definition at line 81 of file LinearIsotropicHardening.hpp.

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

◆ doWork() [2/2]

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

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

Definition at line 91 of file simple_elasticity.cpp.

94  {
95 
97 
98  // get number of dofs on row
99  nbRows = row_data.getIndices().size();
100  // if no dofs on row, exit that work, nothing to do here
101  if (!nbRows)
103 
104  // get number of dofs on column
105  nbCols = col_data.getIndices().size();
106  // if no dofs on Columbia, exit nothing to do here
107  if (!nbCols)
109 
110  // K_ij matrix will have 3 times the number of degrees of freedom of the
111  // i-th entity set (nbRows)
112  // and 3 times the number of degrees of freedom of the j-th entity set
113  // (nbCols)
114  K.resize(nbRows, nbCols, false);
115  K.clear();
116 
117  // get number of integration points
118  nbIntegrationPts = getGaussPts().size2();
119  // check if entity block is on matrix diagonal
120  if (row_side == col_side && row_type == col_type) {
121  isDiag = true;
122  } else {
123  isDiag = false;
124  }
125 
126  // integrate local matrix for entity block
127  CHKERR iNtegrate(row_data, col_data);
128 
129  // assemble local matrix
130  CHKERR aSsemble(row_data, col_data);
131 
133  }
MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Assemble local entity block matrix.
bool isDiag
true if this block is on diagonal
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
MatrixDouble K
const VectorInt & getIndices() const
Get global indices of dofs on entity.
int nbRows
number of dofs on rows
#define CHKERR
Inline error check.
Definition: definitions.h:596
int nbIntegrationPts
number of integration points
int nbCols
number if dof on column
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Integrate B^T D B operator.

◆ iNtegrate() [1/2]

MoFEMErrorCode OpK::iNtegrate ( DataForcesAndSourcesCore::EntData row_data,
DataForcesAndSourcesCore::EntData col_data 
)
protected

Definition at line 133 of file LinearIsotropicHardening.hpp.

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 
151 
153  // get element volume
154  double vol = getVolume();
155 
156  // get intergrayion weights
157  auto t_w = getFTensor0IntegrationWeight();
158 
159  // get derivatives of base functions on rows
160  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
161  // iterate over integration points
162  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
163 
164  // calculate scalar weight times element volume
165  const double a = t_w * vol;
166 
167  // iterate over row base functions
168  for (int rr = 0; rr != nbRows / 3; ++rr) {
169 
170  // get sub matrix for the row
171  auto t_m = get_tensor2(K, 3 * rr, 0);
172 
173  // get derivatives of base functions for columns
174  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
175 
176  // iterate column base functions
177  for (int cc = 0; cc != nbCols / 3;++cc) {
178 
179  // integrate block local stiffens matrix
180  t_m(i, k) += a * (MT_1(i, j, k, l) *
181  (t_row_diff_base(j) * t_col_diff_base(l)));
182 
183  // move to next column base function
184  ++t_col_diff_base;
185 
186  // move to next block of local stiffens matrix
187  ++t_m;
188  }
189 
190  // move to next row base function
191  ++t_row_diff_base;
192  }
193 
194  // move to next integration weight
195  ++t_w;
196  }
197 
199  }
CommonData & commonData
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MatrixDouble consistentMatrix
MatrixDouble K
int nbRows
number of dofs on rows
#define TENSOR4_MAT_PTR(MAT)
Definition: definitions.h:646
int nbIntegrationPts
number of integration points
int nbCols
number if dof on column
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

◆ iNtegrate() [2/2]

MoFEMErrorCode OpK::iNtegrate ( DataForcesAndSourcesCore::EntData row_data,
DataForcesAndSourcesCore::EntData col_data 
)
protected

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

Definition at line 148 of file simple_elasticity.cpp.

149  {
151 
152  // get sub-block (3x3) of local stiffens matrix, here represented by second
153  // order tensor
154  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
156  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
157  &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
158  &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
159  };
160 
165 
166  // get element volume
167  double vol = getVolume();
168 
169  // get intergration weights
170  auto t_w = getFTensor0IntegrationWeight();
171 
172  // get derivatives of base functions on rows
173  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
174  // iterate over integration points
175  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
176 
177  // calculate scalar weight times element volume
178  const double a = t_w * vol;
179 
180  // iterate over row base functions
181  for (int rr = 0; rr != nbRows / 3; ++rr) {
182 
183  // get sub matrix for the row
184  auto t_m = get_tensor2(K, 3 * rr, 0);
185 
186  // get derivatives of base functions for columns
187  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
188 
189  // iterate column base functions
190  for (int cc = 0; cc != nbCols / 3;++cc) {
191 
192  // integrate block local stiffens matrix
193  t_m(i, k) +=
194  a * (tD(i, j, k, l) * (t_row_diff_base(j) * t_col_diff_base(l)));
195 
196  // move to next column base function
197  ++t_col_diff_base;
198 
199  // move to next block of local stiffens matrix
200  ++t_m;
201  }
202 
203  // move to next row base function
204  ++t_row_diff_base;
205  }
206 
207  // move to next integration weight
208  ++t_w;
209  }
210 
212  }
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Ddg< double, 3, 3 > tD
MatrixDouble K
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
int nbCols
number if dof on column
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407

Member Data Documentation

◆ commonData

CommonData& OpK::commonData

Definition at line 64 of file LinearIsotropicHardening.hpp.

◆ isDiag

bool OpK::isDiag
protected

true if this block is on diagonal

Definition at line 139 of file simple_elasticity.cpp.

◆ K

MatrixDouble OpK::K

Definition at line 32 of file simple_elasticity.cpp.

◆ nbCols

int OpK::nbCols
protected

number if dof on column

Definition at line 137 of file simple_elasticity.cpp.

◆ nbIntegrationPts

int OpK::nbIntegrationPts
protected

number of integration points

Definition at line 138 of file simple_elasticity.cpp.

◆ nbRows

int OpK::nbRows
protected

number of dofs on rows

Definition at line 136 of file simple_elasticity.cpp.

◆ pOisson

double OpK::pOisson

Definition at line 40 of file simple_elasticity.cpp.

◆ tD

FTensor::Ddg<double, 3, 3> OpK::tD

Definition at line 35 of file simple_elasticity.cpp.

◆ yOung

double OpK::yOung

Definition at line 38 of file simple_elasticity.cpp.


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