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: elasticity_mixed_formulation.cpp, and ElasticityMixedFormulation.hpp.

Definition at line 311 of file ElasticityMixedFormulation.hpp.

Constructor & Destructor Documentation

◆ OpAssembleK()

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

inline

Examples: ElasticityMixedFormulation.hpp.

Definition at line 321 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: ElasticityMixedFormulation.hpp.

Definition at line 327 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_NULL ? 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);
   }