Added a new HDiv Shape called TPZShapeHDivOptimized

Mesh/TPZCompElHDivDuplConnects.cpp

@@ -2,6 +2,7 @@
 #include "TPZCompElHDivDuplConnectsBound.h"
 #include "TPZMaterial.h"
 #include "TPZShapeHDiv.h"
+#include "TPZShapeHDivOptimized.h"
 #include "TPZShapeHDivConstant.h"
 #include "pzlog.h"
 #include "pzconnect.h"
@@ -108,6 +109,26 @@ int TPZCompElHDivDuplConnects<TSHAPE>::NConnectShapeF(int connect, int order)con
             return nshape;
         }
         break;
+    case HDivFamily::EHDivOptimized:
+        {
+            int conCorrect = connect/2;
+            int res = connect % 2;
+            int nshape;
+            if (!fDuplicationActive){
+                return TPZShapeHDivOptimized<TSHAPE>::ComputeNConnectShapeF(connect,order);
+            } else {
+                nshape = TPZShapeHDivOptimized<TSHAPE>::ComputeNConnectShapeF(conCorrect,order);
+            }
+            if (res == 1){ 
+                nshape -= 1;
+            } else {
+                if (connect != 2*TSHAPE::NFacets){
+                    nshape = 1;
+                }
+            }
+            return nshape;   
+        }
+        break;
 
     default:
         return -1;
@@ -145,6 +166,10 @@ void TPZCompElHDivDuplConnects<TSHAPE>::InitMaterialData(TPZMaterialData &data)
         TPZShapeHDivConstant<TSHAPE>::Initialize(ids, orders, sideorient, data);
         nvec_shape = this->NShapeF();
         break;
+    case HDivFamily::EHDivOptimized:
+        TPZShapeHDivOptimized<TSHAPE>::Initialize(ids, orders, sideorient, data);
+        nvec_shape = TPZShapeHDivOptimized<TSHAPE>::NShapeF(shapedata);
+        break;
 
     default:
         break;

Mesh/pzelchdiv.cpp

@@ -21,6 +21,7 @@
 #include "tpzline.h"
 #include "tpztriangle.h"
 #include "TPZShapeHDiv.h"
+#include "TPZShapeHDivOptimized.h"
 #include "TPZShapeH1.h"
 #include "TPZShapeHDivConstant.h"
 #include "TPZShapeHCurlNoGrads.h"
@@ -255,7 +256,9 @@ int TPZCompElHDiv<TSHAPE>::NConnectShapeF(int connect, int order)const
     case HDivFamily::EHDivConstant:
         return TPZShapeHDivConstant<TSHAPE>::ComputeNConnectShapeF(connect,order);
         break;
-
+    case HDivFamily::EHDivOptimized:
+        return TPZShapeHDivOptimized<TSHAPE>::ComputeNConnectShapeF(connect,order);
+        break;
     default:
         return -1;
         break;
@@ -879,7 +882,10 @@ void TPZCompElHDiv<TSHAPE>::InitMaterialData(TPZMaterialData &data)
         TPZShapeHDivConstant<TSHAPE>::Initialize(ids, orders, sideorient, data);
         nvec_shape = this->NShapeF();
         break;
-
+    case HDivFamily::EHDivOptimized:
+        TPZShapeHDivOptimized<TSHAPE>::Initialize(ids, orders, sideorient, data);
+        nvec_shape = TPZShapeHDivOptimized<TSHAPE>::NShapeF(shapedata);
+        break;
     default:
         break;
     }
@@ -940,13 +946,14 @@ void TPZCompElHDiv<TSHAPE>::ComputeShape(TPZVec<REAL> &qsi, TPZMaterialData &dat
     case HDivFamily::EHDivStandard:
         TPZShapeHDiv<TSHAPE>::Shape(qsi, shapedata, phiMaster, data.divphi);
         break;
-
     case HDivFamily::EHDivConstant:
         phiMaster.Resize(TSHAPE::Dimension,nshape);
         data.divphi.Resize(nshape,1);
         TPZShapeHDivConstant<TSHAPE>::Shape(qsi, shapedata, phiMaster, data.divphi);
         break;
-
+    case HDivFamily::EHDivOptimized:
+        TPZShapeHDivOptimized<TSHAPE>::Shape(qsi, shapedata, phiMaster, data.divphi);
+        break;
     default:
         DebugStop();
         break;
@@ -992,14 +999,16 @@ void TPZCompElHDiv<TSHAPE>::ComputeShape(TPZVec<REAL> &qsi, TPZMaterialData &dat
         switch (fhdivfam)
         {
         case HDivFamily::EHDivStandard:
-                TPZShapeHDiv<TSHAPE>::Shape(qsifad, shapedata, phiMasterFad, divphiFad);
+            TPZShapeHDiv<TSHAPE>::Shape(qsifad, shapedata, phiMasterFad, divphiFad);
             break;
-
         case HDivFamily::EHDivConstant:
             phiMasterFad.Resize(TSHAPE::Dimension,nshape);
             divphiFad.Resize(nshape,1);
             TPZShapeHDivConstant<TSHAPE>::Shape(qsifad, shapedata, phiMasterFad, divphiFad);
             break;
+        case HDivFamily::EHDivOptimized:
+            TPZShapeHDivOptimized<TSHAPE>::Shape(qsifad, shapedata, phiMasterFad, divphiFad);
+            break;
 
         default:
             DebugStop();

Mesh/pzelchdivbound2.cpp

@@ -70,7 +70,7 @@ TPZIntelGen<TSHAPE>(mesh,gel,1), fSideOrient(1), fhdivfam(hdivfam){
 	this->fIntRule.SetOrder(order);
 
 
-    if (fhdivfam == HDivFamily::EHDivConstant) {
+    if (fhdivfam == HDivFamily::EHDivConstant || fhdivfam == HDivFamily::EHDivOptimized) {
         // For HDiv constant, polynomial order was compatibilized in connectorders, 
         // see TPZShapeHDivConstantBound<TSHAPE>::Initialize. So now we need to update
         // the number of shape functions and also the integration rule
@@ -407,6 +407,7 @@ void TPZCompElHDivBound2<TSHAPE>::InitMaterialData(TPZMaterialData &data)
         TPZShapeHDivBound<TSHAPE>::Initialize(id, connectorder, sideorient, data);
         break;
     case HDivFamily::EHDivConstant:
+    case HDivFamily::EHDivOptimized:
         TPZShapeHDivConstantBound<TSHAPE>::Initialize(id, connectorder, sideorient, data);
         break;
 
@@ -580,6 +581,7 @@ void TPZCompElHDivBound2<TSHAPE>::ComputeShape(TPZVec<REAL> &intpoint, TPZMateri
         }
         break;
     case HDivFamily::EHDivConstant:
+    case HDivFamily::EHDivOptimized:
         {
             data.phi.Resize(this->NShapeF(), 1);
             TPZShapeHDivConstantBound<TSHAPE>::Shape(intpoint, shapedata, data.phi);

Pre/TPZHDivApproxCreator.cpp

@@ -77,8 +77,8 @@ void TPZHDivApproxCreator::CheckSetupConsistency() {
         DebugStop();
     }
 
-    if(fHybridType == HybridizationType::ESemi && fHDivFam != HDivFamily::EHDivConstant){
-        std::cout << "The only HDiv space with available Semi hybridization is HDivConstant" << std::endl;
+    if(fHybridType == HybridizationType::ESemi && (fHDivFam != HDivFamily::EHDivConstant && fHDivFam != HDivFamily::EHDivOptimized)){
+        std::cout << "The only HDiv spaces with available Semi hybridization is HDivConstant and HDivOptimized" << std::endl;
         DebugStop();
     }
 

Pre/pzcreateapproxspace.cpp

@@ -537,7 +537,7 @@ void TPZCreateApproximationSpace::SetAllCreateFunctionsHDivDuplConnects(int dime
     fStyle = EHDiv;
     const HDivFamily &hdivfam = this->fhdivfam;
 
-    if (hdivfam != HDivFamily::EHDivConstant){
+    if (hdivfam != HDivFamily::EHDivConstant || hdivfam != HDivFamily::EHDivOptimized){
         std::cout << "HDiv Dupl Connects not implemented yet for this HDiv family!" << std::endl;
         DebugStop();
     }

Shape/CMakeLists.txt

@@ -31,6 +31,7 @@ set(headers
     TPZShapeHDivBound.h
     TPZShapeHDivCollapsed.h
     TPZEnumApproxFamily.h
+    TPZShapeHDivOptimized.h
     )
 
 set(sources
@@ -55,6 +56,7 @@ set(sources
     TPZShapeHCurlNoGrads.cpp
     TPZShapeHDivBound.cpp
     TPZShapeHDivCollapsed.cpp
+    TPZShapeHDivOptimized.cpp
     )
 
 install(FILES ${headers} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/Shape)

Shape/TPZEnumApproxFamily.h

@@ -2,7 +2,7 @@
 #define TPZENUMAPPROXFAMILY_H
 
 /// Enum stating which flavor of HDiv spaces is being used
-enum class HDivFamily {EHDivStandard,EHDivConstant,EHDivKernel};
+enum class HDivFamily {EHDivStandard,EHDivConstant,EHDivKernel,EHDivOptimized};
 
 /// Enum stating which flavor of H1 spaces is being used
 enum class H1Family {EH1Standard};

Shape/TPZShapeData.h

@@ -46,7 +46,7 @@ class TPZShapeData : virtual public TPZSavable {
     /** @brief Prints the data in a format suitable for Mathematica */
     void PrintMathematica(std::ostream &out) const;    
 
-    static constexpr int MatDataNumPhi{60};
+    static constexpr int MatDataNumPhi{200};
     static constexpr int MatDataNumDir{81};
     static constexpr int MatDataDimSol{10};//TODO:Remove?
     static constexpr int MatDataNumSol{20};//TODO:Remove?

Shape/TPZShapeHDiv.cpp

@@ -178,7 +178,7 @@ void TPZShapeHDiv<TSHAPE>::ComputeVecandShape(TPZShapeData &data) {
     // VectorSide indicates the side associated with each vector entry
     TPZManVector<int64_t,27> FirstIndex(TSHAPE::NSides+1);
     // the first index of the shape functions
-    FirstShapeIndex(FirstIndex,scalarorder);
+    FirstShapeIndex(FirstIndex,data.fH1.fConnectOrders);
 
     int64_t nvec = VectorSides.NElements();
     count = 0;
@@ -378,14 +378,14 @@ void TPZShapeHDiv<TSHAPE>::Shape(const TPZVec<Fad<REAL>> &pt, TPZShapeData &data
 
 
 template<class TSHAPE>
-void TPZShapeHDiv<TSHAPE>::FirstShapeIndex(TPZVec<int64_t> &Index, int &scalarorders) {
+void TPZShapeHDiv<TSHAPE>::FirstShapeIndex(TPZVec<int64_t> &Index, const TPZVec<int> &scalarorders) {
     Index[0] = 0;
 
     for(int iside=0;iside<TSHAPE::NSides;iside++)
     {
         int sideorder = 1;
         if (iside >= TSHAPE::NCornerNodes) {
-            sideorder = scalarorders;
+            sideorder = scalarorders[iside-TSHAPE::NCornerNodes];
         }
         int temp = Index[iside] + TSHAPE::NConnectShapeF(iside,sideorder);
         Index[iside+1] = temp;

Shape/TPZShapeHDiv.h

@@ -32,7 +32,7 @@ struct TPZShapeHDiv
 
     static void IndexShapeToVec(TPZShapeData &data);
 
-    static void FirstShapeIndex(TPZVec<int64_t> &Index, int &scalarorders);
+    static void FirstShapeIndex(TPZVec<int64_t> &Index, const TPZVec<int> &scalarorders);
 
     static void FillOrderScalarShapeFunctions(const TPZVec<int> &connectorders, TPZVec<int> &scalarOrder);