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

#include <users_modules/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, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::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,
		DataForcesAndSourcesCore::EntData &	row_data,
		DataForcesAndSourcesCore::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 = getNormalsAtGaussPtsF3(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 = getNormalsAtGaussPtsF3(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);
     }