Inheritance diagram for OpK:

Collaboration diagram for OpK:

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, nonlinear_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.

      : VolumeElementForcesAndSourcesCore::UserDataOperator("U", "U", OPROWCOL,
                                                            symm) {
 
    // Evaluation of the elastic stiffness tensor, D
 
    // hardcoded choice of elastic parameters
    pOisson = 0.1;
    yOung = 10;
 
    // coefficient used in intermediate calculation
    const double coefficient = yOung / ((1 + pOisson) * (1 - 2 * pOisson));
 
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    FTensor::Index<'k', 3> k;
    FTensor::Index<'l', 3> l;
 
    tD(i, j, k, l) = 0.;
 
    tD(0, 0, 0, 0) = 1 - pOisson;
    tD(1, 1, 1, 1) = 1 - pOisson;
    tD(2, 2, 2, 2) = 1 - pOisson;
 
    tD(0, 1, 0, 1) = 0.5 * (1 - 2 * pOisson);
    tD(0, 2, 0, 2) = 0.5 * (1 - 2 * pOisson);
    tD(1, 2, 1, 2) = 0.5 * (1 - 2 * pOisson);
 
    tD(0, 0, 1, 1) = pOisson;
    tD(1, 1, 0, 0) = pOisson;
    tD(0, 0, 2, 2) = pOisson;
    tD(2, 2, 0, 0) = pOisson;
    tD(1, 1, 2, 2) = pOisson;
    tD(2, 2, 1, 1) = pOisson;
 
    tD(i, j, k, l) *= coefficient;
  }

Member Function Documentation

◆ aSsemble()

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

inlineprotected

Assemble local entity block matrix.

Parameters

row_data	row data (consist base functions on row entity)
col_data	column data (consist base functions on column entity)

Returns: error code

Examples: simple_elasticity.cpp.

Definition at line 205 of file simple_elasticity.cpp.

                                                              {
    MoFEMFunctionBegin;
    // get pointer to first global index on row
    const int *row_indices = &*row_data.getIndices().data().begin();
    // get pointer to first global index on column
    const int *col_indices = &*col_data.getIndices().data().begin();
    Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
                                               : getFEMethod()->snes_B;
    // assemble local matrix
    CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
                        &*K.data().begin(), ADD_VALUES);
 
    if (!isDiag && sYmm) {
      // if not diagonal term and since global matrix is symmetric assemble
      // transpose term.
      K = trans(K);
      CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
                          &*K.data().begin(), ADD_VALUES);
    }
    MoFEMFunctionReturn(0);
  }

◆ 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_side	row side number (local number) of entity on element
col_side	column side number (local number) of entity on element
row_type	type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
col_type	type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
row_data	data for row
col_data	data for column

Returns: error code

Examples: simple_elasticity.cpp.

Definition at line 77 of file simple_elasticity.cpp.

                                                            {
 
    MoFEMFunctionBegin;
 
    // get number of dofs on row
    nbRows = row_data.getIndices().size();
    // if no dofs on row, exit that work, nothing to do here
    if (!nbRows)
      MoFEMFunctionReturnHot(0);
 
    // get number of dofs on column
    nbCols = col_data.getIndices().size();
    // if no dofs on Columbia, exit nothing to do here
    if (!nbCols)
      MoFEMFunctionReturnHot(0);
 
    // K_ij matrix will have 3 times the number of degrees of freedom of the
    // i-th entity set (nbRows)
    // and 3 times the number of degrees of freedom of the j-th entity set
    // (nbCols)
    K.resize(nbRows, nbCols, false);
    K.clear();
 
    // get number of integration points
    nbIntegrationPts = getGaussPts().size2();
    // check if entity block is on matrix diagonal
    if (row_side == col_side && row_type == col_type) {
      isDiag = true;
    } else {
      isDiag = false;
    }
 
    // integrate local matrix for entity block
    CHKERR iNtegrate(row_data, col_data);
 
    // assemble local matrix
    CHKERR aSsemble(row_data, col_data);
 
    MoFEMFunctionReturn(0);
  }

◆ iNtegrate()

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

inlineprotected

Integrate B^T D B operator.

Parameters

row_data	row data (consist base functions on row entity)
col_data	column data (consist base functions on column entity)

Returns: error code

Examples: simple_elasticity.cpp.

Definition at line 133 of file simple_elasticity.cpp.

                                                               {
    MoFEMFunctionBegin;
 
    // get sub-block (3x3) of local stiffens matrix, here represented by second
    // order tensor
    auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
      return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
          &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
          &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
          &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
    };
 
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    FTensor::Index<'k', 3> k;
    FTensor::Index<'l', 3> l;
 
    // get element volume
    double vol = getVolume();
 
    // get intergration weights
    auto t_w = getFTensor0IntegrationWeight();
 
    // get derivatives of base functions on rows
    auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
    // iterate over integration points
    for (int gg = 0; gg != nbIntegrationPts; ++gg) {
 
      // calculate scalar weight times element volume
      const double a = t_w * vol;
 
      // iterate over row base functions
      for (int rr = 0; rr != nbRows / 3; ++rr) {
 
        // get sub matrix for the row
        auto t_m = get_tensor2(K, 3 * rr, 0);
 
        // get derivatives of base functions for columns
        auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
 
        // iterate column base functions
        for (int cc = 0; cc != nbCols / 3; ++cc) {
 
          // integrate block local stiffens matrix
          t_m(i, k) +=
              a * (tD(i, j, k, l) * (t_row_diff_base(j) * t_col_diff_base(l)));
 
          // move to next column base function
          ++t_col_diff_base;
 
          // move to next block of local stiffens matrix
          ++t_m;
        }
 
        // move to next row base function
        ++t_row_diff_base;
      }
 
      // move to next integration weight
      ++t_w;
    }
 
    MoFEMFunctionReturn(0);
  }