\biref operator to calculate and assemble Cmat More...

#include "users_modules/mofem_um_homogenisation/src/BCs_RVELagrange_Periodic.hpp"

Inheritance diagram for BCs_RVELagrange_Periodic::OpRVEBCsPeriodicCalAssemCmat:

Collaboration diagram for BCs_RVELagrange_Periodic::OpRVEBCsPeriodicCalAssemCmat:

Public Member Functions
	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)

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

Public Attributes
Mat	Aij

RVEBC_Data_Periodic &	dAta

CommonFunctionsPeriodic	commonFunctionsPeriodic

CommonDataPeriodic &	commonDataPeriodic

bool	hoGeometry

MatrixDouble	Kmat

MatrixDouble	transKmat

ublas::vector< int >	rrvec

ublas::vector< int >	ccvec

MatrixDouble	NTN

Detailed Description

\biref operator to calculate and assemble Cmat

Definition at line 338 of file BCs_RVELagrange_Periodic.hpp.

Constructor & Destructor Documentation

◆ OpRVEBCsPeriodicCalAssemCmat()

BCs_RVELagrange_Periodic::OpRVEBCsPeriodicCalAssemCmat::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`
	)

inline

Definition at line 345 of file BCs_RVELagrange_Periodic.hpp.

     :
 
    FlatPrismElementForcesAndSurcesCore::UserDataOperator(
      lagrang_field_name, field_name, UserDataOperator::OPROWCOL
    ),Aij(aij),
    dAta(data),
    commonFunctionsPeriodic(common_functions_periodic),
    commonDataPeriodic(common_data_periodic),
    hoGeometry(ho_geometry){
      sYmm = false;  //This will make sure to loop over all intities (e.g. for order=2 it will make doWork to loop 16 time)
    }

Member Function Documentation

◆ doWork()

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

inline

Definition at line 369 of file BCs_RVELagrange_Periodic.hpp.

      {
      PetscFunctionBegin;
 
      try {
//        cout<<"Hello from  OpRVEBCsPeriodicLhsAssemCmat "<<endl;
        if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
        if(col_data.getIndices().size()==0) PetscFunctionReturn(0);
        if(row_type == MBEDGE && row_side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
        if(col_type == MBEDGE && col_side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
        if(row_type == MBTRI && row_side == 4) PetscFunctionReturn(0);  //ignore second triangel face (cananical number 4)
        if(col_type == MBTRI && col_side == 4) PetscFunctionReturn(0);  //ignore second triangel face (cananical number 4)
 
 
        //Cmat calculation
        
        int rank = col_data.getFieldDofs()[0]->getNbOfCoeffs();
 
        MatrixDouble N_mat_row;
        MatrixDouble N_mat_col;
 
 
        //1st face
        for(unsigned int gg = 0;gg<col_data.getN().size1();gg++) {
          double area;
          VectorDouble normal_f3 = getNormalsAtGaussPtF3(gg);
          area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
          double val = getGaussPts()(2,gg)*area;
 
          if(row_type == MBVERTEX) {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.RowN[0][row_type][row_side], N_mat_row, 2); CHKERRQ(ierr);
//            cerr << "N_mat_row  "<< N_mat_row <<  endl;
            }
          else {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.RowN[0][row_type][row_side], N_mat_row); CHKERRQ(ierr);}
 
          if(col_type == MBVERTEX) {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.ColN[0][col_type][col_side], N_mat_col, 2); CHKERRQ(ierr);
//            cerr << "N_mat_row  "<< N_mat_row <<  endl;
            }
          else {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.ColN[0][col_type][col_side], N_mat_col); CHKERRQ(ierr);}
          if(gg==0){
//            cerr << "row_type  "<< row_type <<" row_side "<<row_side << "    "<<  endl;
//            cerr << "col_type  "<< col_type <<" col_side "<< col_side<< "    "<<  endl;
            NTN=prod(trans(N_mat_row),N_mat_col);  //we don't need to define size of NTN
            Kmat=-1.0*val*NTN; //we don't need to defien size of K
          }else{
            NTN=prod(trans(N_mat_row),N_mat_col);
            Kmat+=-1.0*val*NTN;
          }
        }
        rrvec=commonDataPeriodic.RowInd[0][row_type][row_side];
        ccvec=commonDataPeriodic.ColInd[0][col_type][col_side];
        ierr = MatSetValues(Aij,rrvec.size(),&rrvec[0],ccvec.size(),&ccvec[0],&Kmat(0,0),ADD_VALUES); CHKERRQ(ierr);
        transKmat = trans(Kmat);
        ierr = MatSetValues(Aij,ccvec.size(),&ccvec[0],rrvec.size(),&rrvec[0],&transKmat(0,0),ADD_VALUES); CHKERRQ(ierr);
 
//        cerr << "rrvec  "<< rrvec <<  endl;
//        cerr << "ccvec  "<< ccvec <<  endl;
//        cerr << "Kmat  "<< Kmat <<  endl;
 
 
        //2nd face
        for(unsigned int gg = 0;gg<col_data.getN().size1();gg++) {
          double area;
          VectorDouble normal_f3 = getNormalsAtGaussPtF3(gg);
          area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
          double val = getGaussPts()(2,gg)*area;
 
          if(row_type == MBVERTEX) {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.RowN[1][row_type][row_side], N_mat_row, 2); CHKERRQ(ierr);
//            cerr << "N_mat_row  "<< N_mat_row <<  endl;
            }
          else {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.RowN[1][row_type][row_side], N_mat_row); CHKERRQ(ierr);}
 
          if(col_type == MBVERTEX) {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.ColN[1][col_type][col_side], N_mat_col, 2); CHKERRQ(ierr);
//            cerr << "N_mat_row  "<< N_mat_row <<  endl;
            }
          else {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.ColN[1][col_type][col_side], N_mat_col); CHKERRQ(ierr);}
          if(gg==0){
//            cerr << "row_type  "<< row_type <<" row_side "<<row_side << "    "<<  endl;
//            cerr << "col_type  "<< col_type <<" col_side "<< col_side<< "    "<<  endl;
            NTN=prod(trans(N_mat_row),N_mat_col);  //we don't need to define size of NTN
            Kmat=val*NTN; //we don't need to defien size of K
          }else{
            NTN=prod(trans(N_mat_row),N_mat_col);
            Kmat+=val*NTN;
          }
        }
        rrvec=commonDataPeriodic.RowInd[1][row_type][row_side];
        ccvec=commonDataPeriodic.ColInd[1][col_type][col_side];
        ierr = MatSetValues(Aij,rrvec.size(),&rrvec[0],ccvec.size(),&ccvec[0],&Kmat(0,0),ADD_VALUES); CHKERRQ(ierr);
        transKmat = trans(Kmat);
        ierr = MatSetValues(Aij,ccvec.size(),&ccvec[0],rrvec.size(),&rrvec[0],&transKmat(0,0),ADD_VALUES); CHKERRQ(ierr);
 
 
//        cerr << "\n\n\nrrvec  "<< rrvec <<  endl;
//        cerr << "ccvec  "<< ccvec <<  endl;
//        cerr << "Kmat  "<< Kmat <<  endl;
 
 
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
//      string wait;
//      cin>>wait;
      PetscFunctionReturn(0);
    }