16 #ifndef __BCS_RVELAGRANGE_PERIIODIC_HPP
17 #define __BCS_RVELAGRANGE_PERIIODIC_HPP
53 const string element_name,
54 const string field_name,
55 const string lagrang_field_name,
56 const string mesh_nodals_positions,
57 Range &periodic_prisms
85 setOfRVEBC[1].tRis = tris.subset_by_type(MBTRI);
87 PetscFunctionReturn(0);
97 int shape_size=col_data.getN().size2()/div;
100 N_mat.resize(rank,shape_size*rank); N_mat.clear();
102 for(
int ii=0; ii<shape_size; ii++){
103 for(
int jj=0; jj<rank; jj++){
105 N_mat(jj,kk)=col_data.getN()(gg,gg1); kk++;
109 PetscFunctionReturn(0);
120 int shape_size=RowN.size2()/div;
123 N_mat.resize(rank,shape_size*rank); N_mat.clear();
125 for(
int ii=0; ii<shape_size; ii++){
126 for(
int jj=0; jj<rank; jj++){
128 N_mat(jj,kk)=RowN(gg,gg1); kk++;
132 PetscFunctionReturn(0);
141 ublas::vector<map<EntityType, map<int, map<EntityType, map<int, MatrixDouble > > > > >
Cmat;
142 ublas::vector<map<EntityType, map<int, ublas::vector<int> > > >
RowInd;
143 ublas::vector<map<EntityType, map<int, ublas::vector<int> > > >
ColInd;
147 ublas::vector<map<EntityType, map<int, MatrixDouble > > >
RowN;
148 ublas::vector<map<EntityType, map<int, MatrixDouble > > >
ColN;
169 const string field_name,
173 bool ho_geometry =
false
198 if(data.getIndices().size()==0) PetscFunctionReturn(0);
202 nb=data.getIndices().size()/2;
206 for(
int ii=0; ii<nb; ii++){
240 SETERRQ(PETSC_COMM_SELF,1,
"data inconsitency");
242 PetscFunctionReturn(0);
257 const string field_name,
261 bool ho_geometry =
false
280 if(data.getIndices().size()==0) PetscFunctionReturn(0);
281 if(
type == MBEDGE && side >= 3) PetscFunctionReturn(0);
282 if(
type == MBTRI && side == 4) PetscFunctionReturn(0);
294 nb=data.getIndices().size()/2;
299 for(
int ii=0; ii<nb; ii++){
327 SETERRQ(PETSC_COMM_SELF,1,
"data inconsitency");
332 PetscFunctionReturn(0);
346 const string field_name,
347 const string lagrang_field_name,
352 bool ho_geometry =
false
370 int row_side,
int col_side,
371 EntityType row_type,EntityType col_type,
379 if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
380 if(col_data.getIndices().size()==0) PetscFunctionReturn(0);
381 if(row_type == MBEDGE && row_side >= 3) PetscFunctionReturn(0);
382 if(col_type == MBEDGE && col_side >= 3) PetscFunctionReturn(0);
383 if(row_type == MBTRI && row_side == 4) PetscFunctionReturn(0);
384 if(col_type == MBTRI && col_side == 4) PetscFunctionReturn(0);
389 int rank = col_data.getFieldDofs()[0]->getNbOfCoeffs();
396 for(
unsigned int gg = 0;gg<col_data.getN().size1();gg++) {
399 area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
400 double val = getGaussPts()(2,gg)*area;
402 if(row_type == MBVERTEX) {
409 if(col_type == MBVERTEX) {
418 NTN=prod(trans(N_mat_row),N_mat_col);
421 NTN=prod(trans(N_mat_row),N_mat_col);
437 for(
unsigned int gg = 0;gg<col_data.getN().size1();gg++) {
440 area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
441 double val = getGaussPts()(2,gg)*area;
443 if(row_type == MBVERTEX) {
450 if(col_type == MBVERTEX) {
459 NTN=prod(trans(N_mat_row),N_mat_col);
462 NTN=prod(trans(N_mat_row),N_mat_col);
478 }
catch (
const std::exception& ex) {
480 ss <<
"throw in method: " << ex.what() << endl;
481 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
485 PetscFunctionReturn(0);
501 const string field_name,
502 const string lagrang_field_name,
506 bool ho_geometry =
false
526 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
528 double x,y,z,x1,y1,z1;
529 for(
unsigned int gg = 0;gg<data.getN().size1();gg++) {
532 area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
533 double val = getGaussPts()(2,gg)*area;
536 x = getHOCoordsAtGaussPtsF3()(gg,0);
537 y = getHOCoordsAtGaussPtsF3()(gg,1);
538 z = getHOCoordsAtGaussPtsF3()(gg,2);
540 x1 = getHOCoordsAtGaussPtsF4()(gg,0);
541 y1 = getHOCoordsAtGaussPtsF4()(gg,1);
542 z1 = getHOCoordsAtGaussPtsF4()(gg,2);
552 X_mat[0].resize(rank,6,
false);
X_mat[0].clear();
558 X_mat[1].resize(rank,6,
false);
X_mat[1].clear();
565 X_mat[0].resize(rank,3,
false);
X_mat[0].clear();
568 X_mat[1].resize(rank,3,
false);
X_mat[1].clear();
572 SETERRQ(PETSC_COMM_SELF,1,
"not implemented");
576 if(
type == MBVERTEX) {
593 PetscFunctionReturn(0);
602 const string field_name,
603 const string lagrang_field_name,
608 bool ho_geometry =
false
611 field_name,lagrang_field_name,data,common_data_periodic,common_functions_periodic,ho_geometry
621 if(data.getIndices().size()==0) PetscFunctionReturn(0);
622 if(
type == MBEDGE && side >= 3) PetscFunctionReturn(0);
623 if(
type == MBTRI && side == 4) PetscFunctionReturn(0);
638 ublas::vector<int> rowvec;
639 if(
type == MBVERTEX) {
640 int nb=data.getIndices().size()/2;
642 for(
int ii=0; ii<nb; ii++){
643 rowvec[ii]=data.getIndices()[ii];
646 rowvec=data.getIndices();
651 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
652 int size = (rank == 3) ? 6 : 3;
653 for(
int ii = 0;ii<size;ii++) {
662 PetscFunctionReturn(0);
676 const string field_name,
677 const string lagrang_field_name,
680 boost::ptr_vector<MethodForForceScaling> &methods_op,
684 bool ho_geometry =
false
687 field_name,lagrang_field_name,data,common_data_periodic,common_functions_periodic,ho_geometry
696 if(data.getIndices().size()==0) PetscFunctionReturn(0);
697 if(
type == MBEDGE && side >= 3) PetscFunctionReturn(0);
698 if(
type == MBTRI && side == 4) PetscFunctionReturn(0);
703 ublas::vector<int> rowvec;
704 if(
type == MBVERTEX) {
705 int nb=data.getIndices().size()/2;
707 for(
int ii=0; ii<nb; ii++){
708 rowvec[ii]=data.getIndices()[ii];
711 rowvec=data.getIndices();
716 scaled_given_strain.resize(6);
735 PetscFunctionReturn(0);
748 const string field_name,
752 bool ho_geometry =
false
771 if(data.getIndices().size()==0) PetscFunctionReturn(0);
772 int nb_gauss_pts = data.getN().size1();
773 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
776 if(
type == MBVERTEX) {
782 for(
int gg = 0;gg<nb_gauss_pts;gg++) {
787 int nb=data.getFieldData().size()/2;
791 for(
int ii=0; ii<nb; ii++){
799 for(
unsigned int gg = 0; gg<nb_gauss_pts;gg++) {
800 if(
type == MBVERTEX) {
805 if(
type == MBVERTEX) {
808 else if(
type == MBEDGE && side < 3){
810 else if(
type == MBEDGE && side >= 6){
812 else if(
type == MBTRI && side == 3){
814 else if(
type == MBTRI && side == 4){
818 PetscFunctionReturn(0);
834 const string field_name,
839 bool ho_geometry =
false
860 if(data.getIndices().size()==0) PetscFunctionReturn(0);
861 if(
type == MBQUAD) PetscFunctionReturn(0);
862 if(
type == MBEDGE && side >= 3) PetscFunctionReturn(0);
863 if(
type == MBTRI && side == 4) PetscFunctionReturn(0);
865 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
867 if(
type == MBVERTEX) {
868 int nb=data.getIndices().size()/2;
870 for(
int ii=0; ii<nb; ii++){
871 row_ind[ii]=data.getIndices()[ii];
875 for(
unsigned int gg = 0;gg<data.getN().size1();gg++) {
878 area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
879 double val = getGaussPts()(2,gg)*area;
881 if(
type == MBVERTEX) {
894 if(
type == MBVERTEX) {
901 PetscFunctionReturn(0);
915 const string field_name,
919 bool ho_geometry =
false
940 if(data.getIndices().size()==0) PetscFunctionReturn(0);
941 if(
type == MBEDGE && side >= 3) PetscFunctionReturn(0);
942 if(
type == MBTRI && side == 4) PetscFunctionReturn(0);
944 int nb_gauss_pts = data.getN().size1();
945 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
948 if(
type == MBVERTEX) {
950 for(
int gg = 0;gg<nb_gauss_pts;gg++) {
954 int nb=data.getFieldData().size()/2;
957 for(
int ii=0; ii<nb; ii++){
963 for(
unsigned int gg = 0;gg<data.getN().size1();gg++) {
964 if(
type == MBVERTEX) {
969 if(
type == MBVERTEX) {
975 PetscFunctionReturn(0);
990 const string field_name,
995 bool ho_geometry =
false
1016 if(data.getIndices().size()==0) PetscFunctionReturn(0);
1017 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
1018 if(
type == MBVERTEX) {
1020 int nb=data.getIndices().size()/2;
1023 for(
int ii=0; ii<nb; ii++){
1024 col_ind[0][ii]=data.getIndices()[ii];
1025 col_ind[1][ii]=data.getIndices()[ii+nb];
1031 for(
unsigned int gg = 0;gg<data.getN().size1();gg++) {
1035 area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
1036 double val = getGaussPts()(2,gg)*area;
1038 if(
type == MBVERTEX) {
1045 if(
type == MBVERTEX) {
1048 else if(
type == MBEDGE && side < 3){
1050 else if(
type == MBEDGE && side >= 6){
1052 else if(
type == MBTRI && side == 3){
1054 else if(
type == MBTRI && side == 4){
1057 if(
type == MBVERTEX) {
1060 else if(
type == MBEDGE && side < 3){
1062 else if(
type == MBEDGE && side >= 6){
1064 else if(
type == MBTRI && side == 3){
1066 else if(
type == MBTRI && side == 4){
1073 if(
type == MBVERTEX) {
1076 else if(
type == MBEDGE && side < 3){
1078 else if(
type == MBEDGE && side >= 6){
1080 else if(
type == MBTRI && side == 3){
1082 else if(
type == MBTRI && side == 4){
1084 PetscFunctionReturn(0);
1093 string lagrang_field_name,
1094 string mesh_nodals_positions,
1099 bool ho_geometry =
false;
1170 PetscFunctionReturn(0);
1177 string field_name,
string lagrang_field_name,
string mesh_nodals_positions,Mat aij,vector<Vec> &
f
1181 bool ho_geometry =
false;
1222 PetscFunctionReturn(0);
1232 const string field_name,
1233 const string lagrang_field_name,
1238 bool ho_geometry =
false
1241 field_name,lagrang_field_name,data,common_data_periodic,common_functions_periodic,ho_geometry
1252 if(data.getIndices().size()==0) PetscFunctionReturn(0);
1253 if(
type == MBEDGE && side >= 3) PetscFunctionReturn(0);
1254 if(
type == MBTRI && side == 4) PetscFunctionReturn(0);
1270 if(
type == MBVERTEX) {
1271 int nb=data.getFieldData().size()/2;
1274 for(
int ii=0; ii<nb; ii++){
1291 int Indices6[6]={0, 1, 2, 3, 4, 5};
1292 int Indices3[3]={0, 1, 2};
1294 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
1305 SETERRQ(PETSC_COMM_SELF,1,
"not implemented");
1308 PetscFunctionReturn(0);
1314 string lagrang_field_name,
1315 string mesh_nodals_positions,
1319 bool ho_geometry =
false;
1331 PetscFunctionReturn(0);
ForcesAndSourcesCore::UserDataOperator UserDataOperator
EntitiesFieldData::EntData EntData
#define CHKERRQ_MOAB(a)
check error code of MoAB function
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual bool check_field(const std::string &name) const =0
check if field is in database
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
constexpr auto VecSetValues
constexpr auto MatSetValues
map< int, RVEBC_Data > setOfRVEBC
maps side set id with appropriate FluxData
MoFEM::Interface & mField
ublas::vector< MatrixDouble > D_mat
ublas::vector< map< EntityType, map< int, ublas::vector< int > > > > RowInd
ublas::vector< map< EntityType, map< int, MatrixDouble > > > ColN
ublas::vector< map< EntityType, map< int, MatrixDouble > > > RowN
ublas::vector< VectorDouble > LagMulAtGaussPts
ublas::vector< map< EntityType, map< int, map< EntityType, map< int, MatrixDouble > > > > > Cmat
ublas::vector< ublas::vector< VectorDouble > > DispAtGaussPts
ublas::vector< map< EntityType, map< int, ublas::vector< int > > > > ColInd
PetscErrorCode shapeMatNew(int rank, unsigned int gg, MatrixDouble &RowN, MatrixDouble &N_mat, int div=1)
PetscErrorCode shapeMat(int rank, unsigned int gg, DataForcesAndSourcesCore::EntData &col_data, MatrixDouble &N_mat, int div=1)
MyPrismFE(MoFEM::Interface &_mField)
\biref operator to calculate the RHS of the constrain for the RVE boundary conditions
CommonFunctionsPeriodic & commonFunctionsPeriodic
CommonDataPeriodic & commonDataPeriodic
OpDmatRhs(const string field_name, const string lagrang_field_name, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
RVEBC_Data_Periodic & dAta
PetscErrorCode calculateDmat(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
ublas::vector< MatrixDouble > X_mat
\biref operator to calculate and assemble Cmat
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
ublas::vector< int > rrvec
CommonFunctionsPeriodic commonFunctionsPeriodic
RVEBC_Data_Periodic & dAta
CommonDataPeriodic & commonDataPeriodic
OpRVEBCsPeriodicCalAssemCmat(const string field_name, const string lagrang_field_name, Mat aij, RVEBC_Data_Periodic &data, CommonFunctionsPeriodic &common_functions_periodic, CommonDataPeriodic &common_data_periodic, bool ho_geometry=false)
ublas::vector< int > ccvec
ublas::vector< ublas::vector< int > > col_ind
OpRVEBCsPeriodicCalCTLam(const string field_name, RVEBC_Data_Periodic &data, Vec f, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
CommonFunctionsPeriodic & commonFunctionsPeriodic
RVEBC_Data_Periodic & dAta
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
CommonDataPeriodic & commonDataPeriodic
ublas::vector< int > row_ind
CommonDataPeriodic & commonDataPeriodic
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
RVEBC_Data_Periodic & dAta
CommonFunctionsPeriodic & commonFunctionsPeriodic
OpRVEBCsPeriodicCalCU(const string field_name, RVEBC_Data_Periodic &data, Vec f, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
RVEBC_Data_Periodic & dAta
CommonFunctionsPeriodic & commonFunctionsPeriodic
CommonDataPeriodic & commonDataPeriodic
ublas::vector< VectorDouble > field_data_nodes
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpRVEBCsPeriodicCalDispAtGaussPts(const string field_name, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
OpRVEBCsPeriodicCalLagMulAtGaussPts(const string field_name, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
CommonDataPeriodic & commonDataPeriodic
VectorDouble field_data_nodes
RVEBC_Data_Periodic & dAta
CommonFunctionsPeriodic & commonFunctionsPeriodic
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
CommonDataPeriodic & commonDataPeriodic
OpRVEBCsPeriodicColInd(const string field_name, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
CommonFunctionsPeriodic & commonFunctionsPeriodic
RVEBC_Data_Periodic & dAta
ublas::vector< VectorDouble > field_data_nodes
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
boost::ptr_vector< MethodForForceScaling > & methodsOp
OpRVEBCsPeriodicRhs_givenStrain(const string field_name, const string lagrang_field_name, Vec f, VectorDouble given_strain, boost::ptr_vector< MethodForForceScaling > &methods_op, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
VectorDouble appliedStrain
OpRVEBCsPeriodicRhs(const string field_name, const string lagrang_field_name, vector< Vec > &f, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
CommonDataPeriodic & commonDataPeriodic
ublas::vector< VectorDouble > field_data_nodes
RVEBC_Data_Periodic & dAta
OpRVEBCsPeriodicRowInd(const string field_name, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
CommonFunctionsPeriodic & commonFunctionsPeriodic
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
\biref operator to calculate the RVE homogenised stress
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpRVEHomoStress(const string field_name, const string lagrang_field_name, Vec stress_homo, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
ublas::vector< VectorDouble > Stress_Homo_elem
map< int, RVEBC_Data_Periodic > setOfRVEBCPrisms
maps side set id with appropriate FluxData
MyPrismFE & getLoopFeRVEBCRhs()
MyPrismFE & getLoopFeRVEBCLhs()
boost::ptr_vector< MethodForForceScaling > methodsOp
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
MyPrismFE & getLoopFeRVEBCStress()
CommonDataPeriodic commonDataPeriodic
CommonFunctionsPeriodic commonFunctionsPeriodic
BCs_RVELagrange_Periodic(MoFEM::Interface &m_field)
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions, Range &periodic_prisms)
MyPrismFE feRVEBCRhsResidual
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &f)
PetscErrorCode setRVEBCsOperatorsNonlinear(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &fvec, Vec f, VectorDouble given_strain)
MyPrismFE & getLoopFeRVEBCRhsResidual()
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
virtual moab::Interface & get_moab()=0
Deprecated interface functions.