Getting GFEM to work

Analysis/TPZLinearAnalysis.cpp

@@ -57,7 +57,7 @@ void TPZLinearAnalysis::AssembleT()
   if(!this->fStructMatrix){
     std::cout<<"Setting default struct matrix: ";
     TPZAutoPointer<TPZStructMatrix> str{nullptr};
-#if defined(USING_MKL) || defined(USING_EIGEN)
+#if defined(USING_MKL)// || defined(USING_EIGEN)
     std::cout<<"sparse(non-symmetric)"<<std::endl;
     str = new TPZSpStructMatrix<TVar>(fCompMesh);
 #else

CMakeLists.txt

@@ -387,11 +387,23 @@ if(USING_UMFPACK)
 		  target_compile_definitions(pz PUBLIC USING_UMFPACK)
 endif(USING_UMFPACK)
 
+# option(USING_EIGEN "Whether the EIGEN3 library will be linked in" OFF)
+# if(USING_EIGEN)
+# 		  find_package(Eigen3 CONFIG)
+# 		  target_link_libraries(pz PUBLIC Eigen3::Eigen)
+# 		  target_compile_definitions(pz PUBLIC USING_EIGEN)
+# endif(USING_EIGEN)
+
+
 option(USING_EIGEN "Whether the EIGEN3 library will be linked in" OFF)
 if(USING_EIGEN)
-		  find_package(Eigen3 CONFIG)
-		  target_link_libraries(pz PUBLIC Eigen3::Eigen)
+      include(cmake/EnableEigen.cmake)
+      enable_eigen(pz)
+      include_directories(/Users/philippedevloo/GitHub/GFemResearch/neopz_GFEM/build-cmake/_deps/eigen-src)
+		  # find_package(Eigen3 CONFIG)
+		  # target_link_libraries(pz PUBLIC Eigen3::Eigen)
 		  target_compile_definitions(pz PUBLIC USING_EIGEN)
+		  target_compile_definitions(pz PUBLIC PZ_USING_EIGEN)
 endif(USING_EIGEN)
 
 
@@ -401,7 +413,6 @@ endif(USING_EIGEN)
 
 
 
-
 include(cmake/InstallPZConfigFile.cmake)
 
 #############EXTRA TARGETS#############

Material/DarcyFlow/TPZDarcyFlow.cpp

@@ -182,6 +182,11 @@ int TPZDarcyFlow::NSolutionVariables(int var) const {
 
 void TPZDarcyFlow::Solution(const TPZMaterialDataT<STATE> &data, int var, TPZVec<STATE> &solOut) {
 
+    if(data.fShapeType == TPZMaterialData::EEmpty) {
+        solOut.Resize(NSolutionVariables(var));
+        solOut.Fill(0.);
+        return;
+    }
     switch (var) {
         case 1: {
             // Solution/Pressure

Material/Elasticity/TPZElasticity2D.cpp

@@ -160,10 +160,13 @@ void TPZElasticity2D::Contribute(const TPZMaterialDataT<STATE> &data,
 	REAL nu1 = 1. - nu;//(1-nu)
 	REAL nu2 = (1.-2.*nu)/2.;
 	REAL F = E/((1.+nu)*(1.-2.*nu));
+    // std::cout << "nu1F = " << nu1*F << " \nnu2F = " << nu2*F << " \nF = " << F << "\nnuF = "<< nu*F << std::endl;
 
 	for( int in = 0; in < phr; in++ ) {
 		du(0,0) = dphi(0,in)*axes(0,0)+dphi(1,in)*axes(1,0);//dvx
 		du(1,0) = dphi(0,in)*axes(0,1)+dphi(1,in)*axes(1,1);//dvy
+        // du(0,0) = dphi(0,in);
+        // du(1,0) = dphi(1,in);
 
         for (int col = 0; col < efc; col++) 
         {
@@ -173,6 +176,8 @@ void TPZElasticity2D::Contribute(const TPZMaterialDataT<STATE> &data,
 		for( int jn = 0; jn < phr; jn++ ) {
 			du(0,1) = dphi(0,jn)*axes(0,0)+dphi(1,jn)*axes(1,0);//dux
 			du(1,1) = dphi(0,jn)*axes(0,1)+dphi(1,jn)*axes(1,1);//duy
+            // du(0,1) = dphi(0,jn);
+            // du(1,1) = dphi(1,jn);
 
 
 			if (fPlaneStress != 1){
@@ -235,8 +240,7 @@ void TPZElasticity2D::ContributeBC(const TPZMaterialDataT<STATE> &data,
 
 	int phr = phi.Rows();
 	short in,jn;
-    auto  bcLoadCases =
-        dynamic_cast<TPZMatLoadCasesBC<STATE>&>(bc);
+    auto  bcLoadCases = dynamic_cast<TPZMatLoadCasesBC<STATE>&>(bc);
     if (ef.Cols() != bcLoadCases.NumLoadCases()) {
         DebugStop();
     }
@@ -274,8 +278,8 @@ void TPZElasticity2D::ContributeBC(const TPZMaterialDataT<STATE> &data,
             for(in = 0 ; in < phr; in++) {
                 for (int il = 0; il<NumLoadCases(); il++)
                 {
-                    ef(2*in,il)   += BIGNUMBER * v2[0] * phi(in,0) * weight;        // forced v2 displacement
-                    ef(2*in+1,il) += BIGNUMBER * v2[1] * phi(in,0) * weight;        // forced v2 displacement
+                    ef(2*in,il)   += BIGNUMBER * v2[2*il+0] * phi(in,0) * weight;        // forced v2 displacement
+                    ef(2*in+1,il) += BIGNUMBER * v2[2*il+1] * phi(in,0) * weight;        // forced v2 displacement
                 }
                 for (jn = 0 ; jn < phi.Rows(); jn++)
                 {
@@ -292,8 +296,8 @@ void TPZElasticity2D::ContributeBC(const TPZMaterialDataT<STATE> &data,
             {
                 for (int il = 0; il <fNumLoadCases; il++) 
                 {
-                    ef(2*in,il) += v2[0] * phi(in,0) * weight;        // force in x direction
-                    ef(2*in+1,il) +=  v2[1] * phi(in,0) * weight;      // force in y direction
+                    ef(2*in,il) += v2[2*il+0] * phi(in,0) * weight;        // force in x direction
+                    ef(2*in+1,il) +=  v2[2*il+1] * phi(in,0) * weight;      // force in y direction
                 }
             }
         }
@@ -306,8 +310,8 @@ void TPZElasticity2D::ContributeBC(const TPZMaterialDataT<STATE> &data,
                 for (int il = 0; il <fNumLoadCases; il++) 
                 {
 //                    const auto &v2 = bcLoadCases.GetBCRhsVal(il);
-                    ef(2*in,il) += v2[0] * phi(in,0) * weight;        // force in x direction
-                    ef(2*in+1,il) += v2[1] * phi(in,0) * weight;      // forced in y direction
+                    ef(2*in,il) += v2[2*il+0] * phi(in,0) * weight;        // force in x direction
+                    ef(2*in+1,il) += v2[2*il+1] * phi(in,0) * weight;      // forced in y direction
                 }
 
                 for (jn = 0 ; jn < phi.Rows(); jn++) {
@@ -407,7 +411,7 @@ void TPZElasticity2D::ContributeBC(const TPZMaterialDataT<STATE> &data,
 
 void TPZElasticity2D::FillDataRequirements(TPZMaterialData &data) const 
 {
-    data.fNeedsSol = true;
+    data.fNeedsSol = false;
     data.fNeedsNormal = false;
 
 }
@@ -819,13 +823,14 @@ void TPZElasticity2D::ContributeVecShape(const TPZMaterialDataT<STATE> &data,
 	efc = ef.Cols();
 	ekr = ek.Rows();
 	ekc = ek.Cols();
+    TPZManVector<STATE,3> force(ff);
 
 	if(fForcingFunction) {            // phi(in, 0) :  node in associated forcing function
 		TPZManVector<STATE> res(3);
 		fForcingFunction(data.x,res);
-		ff[0] = res[0];
-		ff[1] = res[1];
-		ff[2] = res[2];
+		force[0] = res[0];
+		force[1] = res[1];
+		force[2] = res[2];
 	}
 
     REAL E(fE_def), nu(fnu_def);
@@ -859,8 +864,8 @@ void TPZElasticity2D::ContributeVecShape(const TPZMaterialDataT<STATE> &data,
 
         for (int col = 0; col < efc; col++) 
         {
-            ef(in,col) += weight*(   ff[0] * phi(0, in)- dphix_i(0,0)*fPreStressXX - dphix_i(1,0)*fPreStressXY
-                                   + ff[1] * phi(1, in)- dphiy_i(0,0)*fPreStressYY - dphiy_i(1,0)*fPreStressXY);
+            ef(in,col) += weight*(   force[0] * phi(0, in)- dphix_i(0,0)*fPreStressXX - dphix_i(1,0)*fPreStressXY
+                                   + force[1] * phi(1, in)- dphiy_i(0,0)*fPreStressYY - dphiy_i(1,0)*fPreStressXY);
         }		
 		for( int jn = 0; jn < phc; jn++ ) {
 
@@ -901,23 +906,51 @@ void TPZElasticity2D::ContributeVecShapeBC(const TPZMaterialDataT<STATE> &data,
                                            TPZFMatrix<STATE> &ef,
                                            TPZBndCondT<STATE> &bc) {
 
-    const TPZFMatrix<REAL> &phi = data.fH1.fPhi;
+    const TPZFMatrix<REAL> &phi = data.phi;
 
 	const REAL &BIGNUMBER  = TPZMaterial::fBigNumber;
 
 	int phc = phi.Cols();
 	short in,jn;
-	auto  bcLoadCases =
+        auto  bcLoadCases = dynamic_cast<TPZMatLoadCasesBC<STATE>&>(bc);
+    if (ef.Cols() != bcLoadCases.NumLoadCases()) {
+        DebugStop();
+    }
+
+//		In general when the problem is  needed to stablish any convention for ContributeBC implementations
+
+    int nstate = NStateVariables();
+
+    const auto nloads = this->fNumLoadCases;
+    constexpr int nvars = 2;
+    const auto &bcNumLoads =
         dynamic_cast<TPZMatLoadCasesBC<STATE>&>(bc);
+
+    TPZManVector<STATE,10> v2(nvars*nloads);
+    TPZFNMatrix<30,STATE> v1(nvars,1);
+	[&bc = std::as_const(bc),
+     &bcNumLoads = std::as_const(bcNumLoads),
+     &data = std::as_const(data),
+     nvars,nloads]( TPZFMatrix<STATE> &v1, TPZVec<STATE> &v2) {
+        if(bc.HasForcingFunctionBC()){
+            bc.ForcingFunctionBC()(data.x,v2,v1);
+        }else {
+            for(auto l = 0; l < nloads; l++){
+                const auto &val2 = bcNumLoads.GetBCRhsVal(l);
+                for(auto i = 0; i < nvars; i++)
+                    v2[nvars*l+i] = val2[i];
+            }
+            v1 = bc.Val1();
+        }
+    }(v1,v2);
+
+	// auto  bcLoadCases = dynamic_cast<TPZMatLoadCasesBC<STATE>&>(bc);
 	switch (bc.Type()) {
 		case 0 :			// Dirichlet condition
 			for(in = 0 ; in < phc; in++) {
                 for (int il = 0; il <fNumLoadCases; il++) 
                 {
-
-                    const auto &v2 = bcLoadCases.GetBCRhsVal(il);
-
-                    ef(in,il) += weight*BIGNUMBER*(v2[0]*phi(0,in) + v2[1] * phi(1,in));
+                    ef(in,il) += weight*BIGNUMBER*(v2[2*il+0]*phi(0,in) + v2[2*il+1] * phi(1,in));
                 }
 				for (jn = 0 ; jn < phc; jn++) {
 
@@ -931,8 +964,7 @@ void TPZElasticity2D::ContributeVecShapeBC(const TPZMaterialDataT<STATE> &data,
             {
                 for (int il = 0; il <fNumLoadCases; il++) 
                 {
-                    const auto &v2 = bcLoadCases.GetBCRhsVal(il);
-                    ef(in,il)+= weight*(v2[0]*phi(0,in) + v2[1]*phi(1,in));
+                    ef(in,il)+= weight*(v2[2*il+0]*phi(0,in) + v2[2*il+1]*phi(1,in));
                 }
             }
 			break;
@@ -942,8 +974,7 @@ void TPZElasticity2D::ContributeVecShapeBC(const TPZMaterialDataT<STATE> &data,
             {
                 for (int il = 0; il <fNumLoadCases; il++) 
                 {
-                    const auto &v2 = bcLoadCases.GetBCRhsVal(il);
-                    ef(in,il) += weight * (v2[0]*phi(0,in) + v2[1]*phi(1,in));
+                    ef(in,il) += weight * (v2[2*il+0]*phi(0,in) + v2[2*il+1]*phi(1,in));
                 }
 
 				for (jn = 0; jn <phc; jn++) {

Material/Elasticity/TPZElasticity2D.h

@@ -22,6 +22,7 @@ class TPZElasticity2D : public TPZMatBase<STATE,
                                           TPZMatErrorSingleSpace<STATE>,
                                           TPZMatLoadCases<STATE>>
 {
+public:
     using TBase = TPZMatBase<STATE,
                              TPZMatSingleSpaceT<STATE>,
                              TPZMatErrorSingleSpace<STATE>,
@@ -43,6 +44,15 @@ public :
 
     TPZElasticity2D(int id);
 
+    TPZElasticity2D(const TPZElasticity2D &copy) : TBase(copy), fE_def(copy.fE_def), fnu_def(copy.fnu_def),
+    fElasticity(copy.fElasticity), ff(copy.ff), fEover21PlusNu_def(copy.fEover21PlusNu_def),
+    fEover1MinNu2_def(copy.fEover1MinNu2_def), fPreStressXX(copy.fPreStressXX),
+    fPreStressYY(copy.fPreStressYY), fPreStressXY(copy.fPreStressXY), fPreStressZZ(copy.fPreStressZZ),
+    fPlaneStress(copy.fPlaneStress)
+    {
+
+    }
+
     /** @name Elasticity */
     /** @brief Set elasticity parameters */
     void SetElasticity(STATE E, STATE nu)

Material/Elasticity/TPZElasticity3D.cpp

@@ -835,8 +835,17 @@ int TPZElasticity3D::NSolutionVariables(int var) const {
 void TPZElasticity3D::Solution(const TPZMaterialDataT<STATE> &data,
                                int var,TPZVec<STATE> &Solout) {
 	const auto &Sol = data.sol[this->fPostProcIndex];
-    const auto &DSol = data.dsol[this->fPostProcIndex];
+    auto DSol = data.dsol[this->fPostProcIndex];
     const auto &axes = data.axes;
+    if(data.fShapeType == TPZMaterialData::EVecShape){
+        TPZFNMatrix<9,STATE> DSolXY(3,3);
+        for(int i = 0; i<3; i++) {
+            for(int j = 0; j<3; j++) {
+                DSolXY(i,j) = DSol(i*3+j,0);
+            }
+        }
+        DSol = DSolXY;
+    }
     TPZFNMatrix<9,STATE> DSolXY(3,3);
     TPZAxesTools<STATE>::Axes2XYZ(DSol, DSolXY, axes);
 	if(var == TPZElasticity3D::EDisplacement) {

Material/TPZMatLoadCases.cpp

@@ -90,6 +90,11 @@ const TPZVec<TVar>& TPZMatLoadCasesBC<TVar>::GetBCRhsVal(int i) const
     if (valVecSize == 0){
         auto tmp =
             dynamic_cast<const TPZBndCondT<TVar>*>(this);
+        if(!tmp){
+            PZError<<__PRETTY_FUNCTION__;
+            PZError<<"\nERROR: Invalid reference for creating BC.\nAborting...\n";
+            DebugStop();
+        }
         return tmp->Val2();        
     }else{
         return fBCRhsValVec[i];

Matrix/TPZEigenSparseMatrix.cpp

@@ -31,7 +31,7 @@ using namespace std;
 
 template<class TVar>
 TPZEigenSparseMatrix<TVar>::TPZEigenSparseMatrix() : TPZRegisterClassId(&TPZEigenSparseMatrix::ClassId),
-TPZFYsmpMatrix<TVar>()
+TPZSYsmpMatrix<TVar>()
 {
 }
 
@@ -41,7 +41,7 @@ TPZFYsmpMatrix<TVar>()
 
 template<class TVar>
 TPZEigenSparseMatrix<TVar>::TPZEigenSparseMatrix(const int64_t rows,const int64_t cols ) :
-TPZRegisterClassId(&TPZEigenSparseMatrix::ClassId),TPZFYsmpMatrix<TVar>(rows,cols) {
+TPZRegisterClassId(&TPZEigenSparseMatrix::ClassId),TPZSYsmpMatrix<TVar>(rows,cols) {
 #ifdef CONSTRUCTOR
 	cerr << "TPZEigenSparseMatrix(int rows,int cols)\n";
 #endif
@@ -151,7 +151,7 @@ TPZMatrix<TVar> &TPZEigenSparseMatrix<TVar>::operator*=(const TVar val)
 template<class TVar>
 void TPZEigenSparseMatrix<TVar>::Print(const char *title, ostream &out ,const MatrixOutputFormat form) const {
 	// Print the matrix along with a identification title
-    TPZFYsmpMatrix<TVar>::Print(title,out,form);
+    TPZSYsmpMatrix<TVar>::Print(title,out,form);
 }
 
 
@@ -220,6 +220,7 @@ int TPZEigenSparseMatrix<TVar>::Decompose(const DecomposeType dt)
     else {
         DebugStop();
     }
+    return 0;
 }
 /**
  * @brief Solves the linear system using Direct methods
@@ -249,6 +250,7 @@ int TPZEigenSparseMatrix<TVar>::SolveDirect ( TPZFMatrix<TVar>& F , const Decomp
     else{
         DebugStop();
     }
+    return 0;
 }
 
 template<class TVar>
@@ -271,6 +273,7 @@ int TPZEigenSparseMatrix<TVar>::SolveDirect ( TPZFMatrix<TVar>& F , const Decomp
     else{
         DebugStop();
     }
+    return 0;
 }
 
 template<class TVar>
@@ -311,7 +314,7 @@ void TPZEigenSparseMatrix<TVar>::MultAdd(const TPZFMatrix<TVar> &x,const TPZFMat
 
 template<class TVar>
 int TPZEigenSparseMatrix<TVar>::ClassId() const{
-    return Hash("TPZEigenSparseMatrix") ^ TPZFYsmpMatrix<TVar>::ClassId() << 1;
+    return Hash("TPZEigenSparseMatrix") ^ TPZSYsmpMatrix<TVar>::ClassId() << 1;
 }
 
 template class TPZEigenSparseMatrix<double>;

Matrix/TPZEigenSparseMatrix.h

@@ -10,7 +10,7 @@
 #include "pz_config.h"
 #include "pzmatrix.h"
 #include "pzfmatrix.h"
-#include "TPZYSMPMatrix.h"
+#include "TPZSYSMPMatrix.h"
 #include <Eigen/Sparse>
 #include <Eigen/Dense>
 
@@ -27,7 +27,7 @@
  * Defines operations on general sparse matrices stored in the (old) Yale Sparse Matrix Package format.
  */
 template<class TVar>
-class TPZEigenSparseMatrix : public TPZFYsmpMatrix<TVar> {
+class TPZEigenSparseMatrix : public TPZSYsmpMatrix<TVar> {
 
     public :
     typedef Eigen::SparseMatrix<TVar,0,int64_t> SpMat; // declares a column-major sparse matrix type of double

Matrix/pzbasematrix.h

@@ -156,7 +156,11 @@ class TPZBaseMatrix : public TPZSavable{
 
   /** @brief decompose the system of equations acording to the decomposition
    * scheme */
-    virtual int Decompose(const DecomposeType dt) = 0;
+    virtual int Decompose(const DecomposeType dt) {
+      std::cout << "TPZBaseMatrix::Decompose is not implemented\n";
+      DebugStop();
+      return 0;
+    }
 
   /** @brief It prints the matrix data in a MatrixFormat Rows X Cols */
   virtual void Print(const char *name, std::ostream &out = std::cout ,const MatrixOutputFormat form = EFormatted) const = 0;

Matrix/pzmatrix.h

@@ -609,8 +609,8 @@ class TPZMatrix: public TPZBaseMatrix
   virtual int64_t Size() const = 0;
   /** @{ */
   /** @brief Pointer to the beginning of the storage of the matrix*/
-  virtual TVar* &Elem() = 0;
-  virtual const TVar* Elem() const = 0;
+  virtual TVar* &Elem() { DebugStop(); static TVar* t{nullptr}; return t; }
+  virtual const TVar* Elem() const { DebugStop(); return nullptr;}
   /** @} */
 	/**
 	 * @brief Is an auxiliar method used by MultiplyAdd