Assemble K *. More...

#include <users_modules/topology_optimization/src/Header.hpp>

Inheritance diagram for OpAssembleK:

Collaboration diagram for OpAssembleK:

Public Member Functions
	OpAssembleK (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...

	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	rowB

MatrixDouble	colB

MatrixDouble	CB

MatrixDouble	K

MatrixDouble	transK

CommonData &	commonData

MatrixDouble	locK

MatrixDouble	translocK

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

DataAtIntegrationPts &	commonData

BlockData &	dAta

Protected Member Functions
virtual MoFEMErrorCode	iNtegrate (DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
	Integrate grad-grad operator. More...

virtual MoFEMErrorCode	aSsemble (DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
	Assemble local entity block matrix. More...

Protected Attributes
int	nbRows
	< error code 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 K *.

*

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

Examples: ElasticityMixedFormulation.hpp, and elasticity_mixed_formulation.cpp.

Definition at line 338 of file Header.hpp.

Constructor & Destructor Documentation

◆ OpAssembleK() [1/2]

OpAssembleK::OpAssembleK	(	CommonData &	common_data,
		bool	symm = `true`
	)

inline

Definition at line 348 of file Header.hpp.

      : VolumeElementForcesAndSourcesCore::UserDataOperator(
            "DISPLACEMENTS", "DISPLACEMENTS", OPROWCOL, symm),
        commonData(common_data) {}

◆ OpAssembleK() [2/2]

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

inline

Definition at line 321 of file ElasticityMixedFormulation.hpp.

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

Member Function Documentation

◆ aSsemble()

virtual MoFEMErrorCode OpAssembleK::aSsemble	(	DataForcesAndSourcesCore::EntData &	row_data,
		DataForcesAndSourcesCore::EntData &	col_data
	)

inlineprotectedvirtual

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

Definition at line 486 of file Header.hpp.

                                                                             {
    MoFEMFunctionBegin;
 
    int nb_dofs_row = row_data.getFieldData().size();
    if (nb_dofs_row == 0)
      MoFEMFunctionReturnHot(0);
    int nb_dofs_col = col_data.getFieldData().size();
    if (nb_dofs_col == 0)
      MoFEMFunctionReturnHot(0);
 
    // 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;
 
    CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
                        &*K.data().begin(), ADD_VALUES);
    if (!isDiag && sYmm) {
      K = trans(K);
      CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
                          &*K.data().begin(), ADD_VALUES);
    }
    MoFEMFunctionReturn(0);
  }

◆ doWork() [1/2]

MoFEMErrorCode OpAssembleK::doWork	(	int	row_side,
		int	col_side,
		EntityType	row_type,
		EntityType	col_type,
		DataForcesAndSourcesCore::EntData &	row_data,
		DataForcesAndSourcesCore::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: ElasticityMixedFormulation.hpp.

Definition at line 363 of file Header.hpp.

                                                                   {
    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);
    // get number of integration points
    nbIntegrationPts = getGaussPts().size2();
    // chekk if entity block is on matrix diagonal
    if (row_side == col_side && row_type == col_type) { // ASK
      isDiag = true;                                    // yes, it is
    } else {
      isDiag = false;
    }
 
    // integrate local matrix for entity block
    CHKERR iNtegrate(row_data, col_data);
    // asseble local matrix
    CHKERR aSsemble(row_data, col_data);
 
    MoFEMFunctionReturn(0);
  }

◆ doWork() [2/2]

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

inline

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);
  }

◆ iNtegrate()

virtual MoFEMErrorCode OpAssembleK::iNtegrate	(	DataForcesAndSourcesCore::EntData &	row_data,
		DataForcesAndSourcesCore::EntData &	col_data
	)

inlineprotectedvirtual

Integrate grad-grad 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

Definition at line 410 of file Header.hpp.

                                                       {
    MoFEMFunctionBegin;
 
    int nb_dofs_row = row_data.getFieldData().size();
    if (nb_dofs_row == 0)
      MoFEMFunctionReturnHot(0);
    int nb_dofs_col = col_data.getFieldData().size();
    if (nb_dofs_col == 0)
      MoFEMFunctionReturnHot(0);
 
    K.resize(nb_dofs_row, nb_dofs_col, false);
    K.clear();
    auto max = [](double a, double b) { return a > b ? a : b; };
    auto min = [](double a, double b) { return a < b ? a : b; };
    auto rho = getFTensor0FromVec(*commonData.dataRhotGaussPtr);
 
          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 intergrayion weights
      auto t_w = getFTensor0IntegrationWeight();
 
      auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
      const int row_nb_base_fun = row_data.getN().size2();
 
      FTensor::Ddg<FTensor::PackPtr<double *, 1>, 3, 3> t_D(
          MAT_TO_DDG(&commonData.D));
 
      for (int gg = 0; gg != nbIntegrationPts; gg++) {
 
        const double coeff = pow(rho, commonData.penalty);
        double a = t_w * vol;
        if (getHoGaussPtsDetJac().size()) {
          a *= getHoGaussPtsDetJac()[gg];
        }
        a *= coeff;
        int rr = 0;
        for (; rr != nbRows / 3; ++rr) {
          auto t_m = get_tensor2(K, 3 * rr, 0);
          auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
          FTensor::Christof<double, 3, 3> t_rowD;
          t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
          for (int cc = 0; cc != nbCols / 3; ++cc) {
            t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
            ++t_col_diff_base;
            ++t_m;
          }
          ++t_row_diff_base;
        }
        for (; rr != row_nb_base_fun; ++rr)
          ++t_row_diff_base;
        ++t_w;
        ++rho;
      }
      MoFEMFunctionReturn(0);
  }

Member Data Documentation

◆ CB

MatrixDouble OpAssembleK::CB

Definition at line 343 of file Header.hpp.

◆ colB

MatrixDouble OpAssembleK::colB

Definition at line 342 of file Header.hpp.

◆ commonData [1/2]

CommonData& OpAssembleK::commonData

Examples: ElasticityMixedFormulation.hpp.

Definition at line 346 of file Header.hpp.

◆ commonData [2/2]

DataAtIntegrationPts& OpAssembleK::commonData

Definition at line 318 of file ElasticityMixedFormulation.hpp.

◆ dAta

BlockData& OpAssembleK::dAta

Examples: ElasticityMixedFormulation.hpp.

Definition at line 319 of file ElasticityMixedFormulation.hpp.

◆ diffDiff

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

Examples: ElasticityMixedFormulation.hpp.

Definition at line 316 of file ElasticityMixedFormulation.hpp.

◆ isDiag

bool OpAssembleK::isDiag

protected

true if this block is on diagonal

Definition at line 401 of file Header.hpp.

◆ K

MatrixDouble OpAssembleK::K

Definition at line 344 of file Header.hpp.

◆ locK

MatrixDouble OpAssembleK::locK

Examples: ElasticityMixedFormulation.hpp.

Definition at line 314 of file ElasticityMixedFormulation.hpp.

◆ nbCols

int OpAssembleK::nbCols

protected

number if dof on column

Definition at line 399 of file Header.hpp.

◆ nbIntegrationPts

int OpAssembleK::nbIntegrationPts

protected

number of integration points

Definition at line 400 of file Header.hpp.

◆ nbRows

int OpAssembleK::nbRows

protected

< error code

number of dofs on rows

Definition at line 398 of file Header.hpp.

◆ rowB

MatrixDouble OpAssembleK::rowB

Definition at line 341 of file Header.hpp.

◆ transK

MatrixDouble OpAssembleK::transK

Definition at line 344 of file Header.hpp.

◆ translocK

MatrixDouble OpAssembleK::translocK

Examples: ElasticityMixedFormulation.hpp.

Definition at line 315 of file ElasticityMixedFormulation.hpp.

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

users_modules/topology_optimization/src/Header.hpp
users_modules/tutorials/cor-7/src/ElasticityMixedFormulation.hpp

Public Member Functions

Public Attributes

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ OpAssembleK() [1/2]

◆ OpAssembleK() [2/2]

Member Function Documentation

◆ aSsemble()

◆ doWork() [1/2]

◆ doWork() [2/2]

◆ iNtegrate()

Member Data Documentation

◆ CB

◆ colB

◆ commonData [1/2]

◆ commonData [2/2]

◆ dAta

◆ diffDiff

◆ isDiag

◆ K

◆ locK

◆ nbCols

◆ nbIntegrationPts

◆ nbRows

◆ rowB

◆ transK

◆ translocK