add support for high-order elements

Antennas/Microwave.pro

@@ -9,6 +9,10 @@ Function{
   DefineFunction[ epsilon, sigma, nu ];
   DefineConstant[ ZL ];
   DefineConstant[ Flag_3Dmodel, Flag_Axisymmetry, Flag_SilverMuller ];
+  DefineConstant[
+    FEorder = {0.5, Choices{0.5="Lowest",1="1st order",1.5="Reduced 2nd order",2="2nd order"},
+      Name "Input/FE order"}
+  ];
 }
 
 If(Flag_3Dmodel)
@@ -27,7 +31,7 @@ EndIf
 
 If(Flag_Axisymmetry)
   Jacobian{
-    { Name JVol ; Case { { Region All ; Jacobian VolAxiSqu ; } } } // or VolAxi
+    { Name JVol ; Case { { Region All ; Jacobian VolAxi ; } } }
     { Name JSur ; Case { { Region All ; Jacobian SurAxi ; } } }
   }
 Else
@@ -43,12 +47,12 @@ Integration {
       { Type Gauss ;
         Case {
 	  { GeoElement Point       ; NumberOfPoints  1 ; }
-	  { GeoElement Line        ; NumberOfPoints  3 ; }
-	  { GeoElement Triangle    ; NumberOfPoints  4 ; }
-	  { GeoElement Quadrangle  ; NumberOfPoints  4 ; }
-	  { GeoElement Tetrahedron ; NumberOfPoints  4 ; }
-	  { GeoElement Hexahedron  ; NumberOfPoints  6 ; }
-	  { GeoElement Prism       ; NumberOfPoints  6 ; }
+	  { GeoElement Line        ; NumberOfPoints  4 ; }
+	  { GeoElement Triangle    ; NumberOfPoints  7 ; }
+	  { GeoElement Quadrangle  ; NumberOfPoints  7 ; }
+	  { GeoElement Tetrahedron ; NumberOfPoints  15 ; }
+	  { GeoElement Hexahedron  ; NumberOfPoints  34; }
+	  { GeoElement Prism       ; NumberOfPoints  21 ; }
 	}
       }
     }
@@ -78,9 +82,29 @@ FunctionSpace {
     BasisFunction {
       { Name se; NameOfCoef ee; Function BF_Edge;
         Support DomainTot ; Entity EdgesOf[All]; }
+      If(FEorder >= 1)
+        { Name se2; NameOfCoef ee2; Function BF_Edge_2E;
+          Support DomainTot ; Entity EdgesOf[All]; }
+      EndIf
+      If(FEorder >= 1.5)
+        { Name se3; NameOfCoef ee3; Function BF_Edge_3F_a;
+          Support DomainTot ; Entity FacetsOf[All]; }
+        { Name se4; NameOfCoef ee4; Function BF_Edge_3F_b;
+          Support DomainTot ; Entity FacetsOf[All]; }
+      EndIf
+      If(FEorder >= 2)
+        { Name se5; NameOfCoef ee5; Function BF_Edge_4E;
+          Support DomainTot ; Entity EdgesOf[All]; }
+      EndIf
     }
     Constraint {
       { NameOfCoef ee; EntityType EdgesOf ; NameOfConstraint ElectricField; }
+      If(FEorder >= 1)
+        { NameOfCoef ee2; EntityType EdgesOf ; NameOfConstraint ElectricField; }
+      EndIf
+      If(FEorder >= 2)
+        { NameOfCoef ee5; EntityType EdgesOf ; NameOfConstraint ElectricField; }
+      EndIf
     }
   }
 
@@ -252,6 +276,10 @@ Resolution {
     }
     Operation {
       CreateDir[Str[myDir]];
+      If(FEorder >= 2)
+        SetGlobalSolverOptions["-petsc_prealloc 400"];
+      EndIf
+
       Generate A; Solve A; SaveSolution A;
       If(Flag_AnalysisType==0)
         PostOperation[Microwave_e];
@@ -340,6 +368,7 @@ PostOperation {
   { Name Microwave_e ; NameOfPostProcessing Microwave_e ;
     Operation {
       Print[ e,  OnElementsOf Region[{Domain,-Pml}], File StrCat[myDir, Sprintf("e_pml%g.pos", !Flag_SilverMuller)] ] ;
+      //Print[ testdR,  OnElementsOf SkinFeed, File "dR.pos" ] ;
       Print[ h_from_e, OnElementsOf Region[{Domain,-Pml}], File StrCat[myDir,Sprintf("h_pml%g.pos", !Flag_SilverMuller)] ];
       Print[ exh,  OnElementsOf Region[{Domain,-Pml}], File StrCat[myDir,Sprintf("exh_pml%g.pos", !Flag_SilverMuller)] ];