v0.13.2
Loading...
Searching...
No Matches
Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | List of all members
OpAssembleK Struct Reference

Assemble K *. More...

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

Inheritance diagram for OpAssembleK:
[legend]
Collaboration diagram for OpAssembleK:
[legend]

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
 
CommonDatacommonData
 
MatrixDouble locK
 
MatrixDouble translocK
 
FTensor::Tensor2< double, 3, 3 > diffDiff
 
DataAtIntegrationPtscommonData
 
BlockDatadAta
 

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.

349 : VolumeElementForcesAndSourcesCore::UserDataOperator(
350 "DISPLACEMENTS", "DISPLACEMENTS", OPROWCOL, symm),
351 commonData(common_data) {}
CommonData & commonData
Definition: Header.hpp:346

◆ OpAssembleK() [2/2]

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

Definition at line 321 of file ElasticityMixedFormulation.hpp.

323 : VolumeElementForcesAndSourcesCore::UserDataOperator("U", "U", OPROWCOL,
324 symm),
325 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_datarow data (consist base functions on row entity)
col_datacolumn data (consist base functions on column entity)
Returns
error code

Definition at line 486 of file Header.hpp.

487 {
489
490 int nb_dofs_row = row_data.getFieldData().size();
491 if (nb_dofs_row == 0)
493 int nb_dofs_col = col_data.getFieldData().size();
494 if (nb_dofs_col == 0)
496
497 // get pointer to first global index on row
498 const int *row_indices = &*row_data.getIndices().data().begin();
499 // get pointer to first global index on column
500 const int *col_indices = &*col_data.getIndices().data().begin();
501 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
502 : getFEMethod()->snes_B;
503
504 CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
505 &*K.data().begin(), ADD_VALUES);
506 if (!isDiag && sYmm) {
507 K = trans(K);
508 CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
509 &*K.data().begin(), ADD_VALUES);
510 }
512 }
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
int nbCols
number if dof on column
Definition: Header.hpp:399
bool isDiag
true if this block is on diagonal
Definition: Header.hpp:401
MatrixDouble K
Definition: Header.hpp:344
int nbRows
< error code
Definition: Header.hpp:398

◆ 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_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
Examples
ElasticityMixedFormulation.hpp.

Definition at line 363 of file Header.hpp.

366 {
368 // get number of dofs on row
369 nbRows = row_data.getIndices().size();
370 // if no dofs on row, exit that work, nothing to do here
371 if (!nbRows)
373 // get number of dofs on column
374 nbCols = col_data.getIndices().size();
375 // if no dofs on Columbia, exit nothing to do here
376 if (!nbCols)
378 // get number of integration points
379 nbIntegrationPts = getGaussPts().size2();
380 // chekk if entity block is on matrix diagonal
381 if (row_side == col_side && row_type == col_type) { // ASK
382 isDiag = true; // yes, it is
383 } else {
384 isDiag = false;
385 }
386
387 // integrate local matrix for entity block
388 CHKERR iNtegrate(row_data, col_data);
389 // asseble local matrix
390 CHKERR aSsemble(row_data, col_data);
391
393 }
virtual MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Integrate grad-grad operator.
Definition: Header.hpp:410
int nbIntegrationPts
number of integration points
Definition: Header.hpp:400
virtual MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Assemble local entity block matrix.
Definition: Header.hpp:486

◆ 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.

330 {
332
333 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
335 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
336 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
337 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
338 };
339
340 const int row_nb_dofs = row_data.getIndices().size();
341 if (!row_nb_dofs)
343 const int col_nb_dofs = col_data.getIndices().size();
344 if (!col_nb_dofs)
346 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
347 dAta.tEts.end()) {
349 }
350 commonData.getBlockData(dAta);
351
352 const bool diagonal_block =
353 (row_type == col_type) && (row_side == col_side);
354 // get number of integration points
355 // Set size can clear local tangent matrix
356 locK.resize(row_nb_dofs, col_nb_dofs, false);
357 locK.clear();
358
359 const int row_nb_gauss_pts = row_data.getN().size1();
360 const int row_nb_base_functions = row_data.getN().size2();
361
362 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
363
368
369 // integrate local matrix for entity block
370 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
371
372 // Get volume and integration weight
373 double w = getVolume() * getGaussPts()(3, gg);
374
375 int row_bb = 0;
376 for (; row_bb != row_nb_dofs / 3; row_bb++) {
377
378 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
379 const int final_bb = diagonal_block ? row_bb + 1 : col_nb_dofs / 3;
380 int col_bb = 0;
381 for (; col_bb != final_bb; col_bb++) {
382
383 auto t_assemble = get_tensor2(locK, 3 * row_bb, 3 * col_bb);
384
385 diffDiff(j, l) =
386 w * row_diff_base_functions(j) * col_diff_base_functions(l);
387
388 t_assemble(i, k) += diffDiff(j, l) * commonData.tD(i, j, k, l);
389 // Next base function for column
390 ++col_diff_base_functions;
391 }
392
393 ++row_diff_base_functions;
394 }
395 for (; row_bb != row_nb_base_functions; row_bb++) {
396 ++row_diff_base_functions;
397 }
398 }
399
400 if (diagonal_block) {
401 for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
402 int col_bb = 0;
403 for (; col_bb != row_bb + 1; col_bb++) {
404 auto t_assemble = get_tensor2(locK, 3 * row_bb, 3 * col_bb);
405 auto t_off_side = get_tensor2(locK, 3 * col_bb, 3 * row_bb);
406 t_off_side(i, k) = t_assemble(k, i);
407 }
408 }
409 }
410
411 const int *row_ind = &*row_data.getIndices().begin();
412 const int *col_ind = &*col_data.getIndices().begin();
413 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
414 : getFEMethod()->ksp_B;
415 CHKERR MatSetValues(B, row_nb_dofs, row_ind, col_nb_dofs, col_ind,
416 &locK(0, 0), ADD_VALUES);
417
418 if (row_type != col_type || row_side != col_side) {
419 translocK.resize(col_nb_dofs, row_nb_dofs, false);
420 noalias(translocK) = trans(locK);
421 CHKERR MatSetValues(B, col_nb_dofs, col_ind, row_nb_dofs, row_ind,
422 &translocK(0, 0), ADD_VALUES);
423 }
424
426 }
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
double w(const double g, const double t)
int r
Definition: sdf.py:5
FTensor::Tensor2< double, 3, 3 > diffDiff

◆ iNtegrate()

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

Integrate grad-grad 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 410 of file Header.hpp.

411 {
413
414 int nb_dofs_row = row_data.getFieldData().size();
415 if (nb_dofs_row == 0)
417 int nb_dofs_col = col_data.getFieldData().size();
418 if (nb_dofs_col == 0)
420
421 K.resize(nb_dofs_row, nb_dofs_col, false);
422 K.clear();
423 auto max = [](double a, double b) { return a > b ? a : b; };
424 auto min = [](double a, double b) { return a < b ? a : b; };
425 auto rho = getFTensor0FromVec(*commonData.dataRhotGaussPtr);
426
427 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
429 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2),
430 &m(r + 1, c + 0), &m(r + 1, c + 1), &m(r + 1, c + 2),
431 &m(r + 2, c + 0), &m(r + 2, c + 1), &m(r + 2, c + 2));
432 };
433
438
439 // get element volume
440 double vol = getVolume();
441
442 // get intergrayion weights
443 auto t_w = getFTensor0IntegrationWeight();
444
445 auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
446 const int row_nb_base_fun = row_data.getN().size2();
447
450
451 for (int gg = 0; gg != nbIntegrationPts; gg++) {
452
453 const double coeff = pow(rho, commonData.penalty);
454 double a = t_w * vol;
455 if (getHoGaussPtsDetJac().size()) {
456 a *= getHoGaussPtsDetJac()[gg];
457 }
458 a *= coeff;
459 int rr = 0;
460 for (; rr != nbRows / 3; ++rr) {
461 auto t_m = get_tensor2(K, 3 * rr, 0);
462 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
464 t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
465 for (int cc = 0; cc != nbCols / 3; ++cc) {
466 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
467 ++t_col_diff_base;
468 ++t_m;
469 }
470 ++t_row_diff_base;
471 }
472 for (; rr != row_nb_base_fun; ++rr)
473 ++t_row_diff_base;
474 ++t_w;
475 ++rho;
476 }
478 }
#define MAT_TO_DDG(SM)
constexpr double a
double rho
Definition: plastic.cpp:101
boost::shared_ptr< VectorDouble > dataRhotGaussPtr
Definition: Header.hpp:28
int penalty
Definition: Header.hpp:47
MatrixDouble D
Definition: Header.hpp:31

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

◆ diffDiff

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

◆ 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

◆ 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

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