Assemble K *. More...

#include "tutorials/cor-7/src/ElasticityMixedFormulation.hpp"

Inheritance diagram for OpAssembleK:

Collaboration diagram for OpAssembleK:

Public Member Functions
	OpAssembleK (DataAtIntegrationPts &common_data, BlockData &data, bool symm=true)

PetscErrorCode	doWork (int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)

Public Attributes
MatrixDouble	locK

MatrixDouble	translocK

FTensor::Tensor2< double, 3, 3 >	diffDiff

DataAtIntegrationPts &	commonData

BlockData &	dAta

Detailed Description

Assemble K *.

*

\[ {\bf{K}} = \int\limits_\Omega {{{\bf{B}}^T}{{\bf{D}}_d}{\bf{B}}d\Omega } \]

Examples: mofem/tutorials/cor-7/elasticity_mixed_formulation.cpp, and mofem/tutorials/cor-7/src/ElasticityMixedFormulation.hpp.

Definition at line 310 of file ElasticityMixedFormulation.hpp.

Constructor & Destructor Documentation

◆ OpAssembleK()

OpAssembleK::OpAssembleK	(	DataAtIntegrationPts &	common_data,
		BlockData &	data,
		bool	symm = `true`
	)

inline

Definition at line 320 of file ElasticityMixedFormulation.hpp.

      : VolumeElementForcesAndSourcesCore::UserDataOperator("U", "U", OPROWCOL,
                                                            symm),
        commonData(common_data), dAta(data) {}

Member Function Documentation

◆ doWork()

PetscErrorCode OpAssembleK::doWork	(	int	row_side,
		int	col_side,
		EntityType	row_type,
		EntityType	col_type,
		EntitiesFieldData::EntData &	row_data,
		EntitiesFieldData::EntData &	col_data
	)

inline

Examples: mofem/tutorials/cor-7/src/ElasticityMixedFormulation.hpp.

Definition at line 326 of file ElasticityMixedFormulation.hpp.

                                                            {
    MoFEMFunctionBegin;
 
    auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
      return FTensor::Tensor2<double *, 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));
    };
 
    const int row_nb_dofs = row_data.getIndices().size();
    if (!row_nb_dofs)
      MoFEMFunctionReturnHot(0);
    const int col_nb_dofs = col_data.getIndices().size();
    if (!col_nb_dofs)
      MoFEMFunctionReturnHot(0);
    if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
        dAta.tEts.end()) {
      MoFEMFunctionReturnHot(0);
    }
    commonData.getBlockData(dAta);
 
    const bool diagonal_block =
        (row_type == col_type) && (row_side == col_side);
    // get number of integration points
    // Set size can clear local tangent matrix
    locK.resize(row_nb_dofs, col_nb_dofs, false);
    locK.clear();
 
    const int row_nb_gauss_pts = row_data.getN().size1();
    const int row_nb_base_functions = row_data.getN().size2();
 
    auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
 
    FTensor::Index<'i', 3> i;
    FTensor::Index<'j', 3> j;
    FTensor::Index<'k', 3> k;
    FTensor::Index<'l', 3> l;
 
    // integrate local matrix for entity block
    for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
 
      // Get volume and integration weight
      double w = getVolume() * getGaussPts()(3, gg);
 
      int row_bb = 0;
      for (; row_bb != row_nb_dofs / 3; row_bb++) {
 
        auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
        const int final_bb = diagonal_block ? row_bb + 1 : col_nb_dofs / 3;
        int col_bb = 0;
        for (; col_bb != final_bb; col_bb++) {
 
          auto t_assemble = get_tensor2(locK, 3 * row_bb, 3 * col_bb);
 
          diffDiff(j, l) =
              w * row_diff_base_functions(j) * col_diff_base_functions(l);
 
          t_assemble(i, k) += diffDiff(j, l) * commonData.tD(i, j, k, l);
          // Next base function for column
          ++col_diff_base_functions;
        }
 
        ++row_diff_base_functions;
      }
      for (; row_bb != row_nb_base_functions; row_bb++) {
        ++row_diff_base_functions;
      }
    }
 
    if (diagonal_block) {
      for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
        int col_bb = 0;
        for (; col_bb != row_bb + 1; col_bb++) {
          auto t_assemble = get_tensor2(locK, 3 * row_bb, 3 * col_bb);
          auto t_off_side = get_tensor2(locK, 3 * col_bb, 3 * row_bb);
          t_off_side(i, k) = t_assemble(k, i);
        }
      }
    }
 
    const int *row_ind = &*row_data.getIndices().begin();
    const int *col_ind = &*col_data.getIndices().begin();
    Mat B = getFEMethod()->ksp_B != PETSC_NULLPTR ? getFEMethod()->ksp_B
                                               : getFEMethod()->ksp_B;
    CHKERR MatSetValues(B, row_nb_dofs, row_ind, col_nb_dofs, col_ind,
                        &locK(0, 0), ADD_VALUES);
 
    if (row_type != col_type || row_side != col_side) {
      translocK.resize(col_nb_dofs, row_nb_dofs, false);
      noalias(translocK) = trans(locK);
      CHKERR MatSetValues(B, col_nb_dofs, col_ind, row_nb_dofs, row_ind,
                          &translocK(0, 0), ADD_VALUES);
    }
 
    MoFEMFunctionReturn(0);
  }