up

7c8956f7 · Christophe Geuzaine · 6a87b110 · 7c8956f7 · 7c8956f7 · 7c8956f7
Commit 7c8956f7 authored 7 years ago by Christophe Geuzaine
--- a/Elasticity/wrench2D.pro
+++ b/Elasticity/wrench2D.pro
@@ -162,14 +162,14 @@ Group {
 Flag_Degree = DefineNumber[ 0, Name "Geometry/Use degree 2 (hierarch.)",
                            Choices{0,1}, Visible 1];
-FE_Degree = ( Flag_Degree == 0 ) ? 1 : 2; // Convert flag value into polynomial degree
+FE_Order = ( Flag_Degree == 0 ) ? 1 : 2; // Convert flag value into polynomial degree
 FunctionSpace {
  { Name H_ux_Mec ; Type Form0 ;
    BasisFunction {
      { Name sxn ; NameOfCoef uxn ; Function BF_Node ;
        Support Dom_H_u_Mec ; Entity NodesOf[ All ] ; }
-     If ( FE_Degree == 2 )
+     If ( FE_Order == 2 )
        { Name sxn2 ; NameOfCoef uxn2 ; Function BF_Node_2E ;
          Support Dom_H_u_Mec; Entity EdgesOf[ All ] ; }
     EndIf
@@ -177,7 +177,7 @@ FunctionSpace {
    Constraint {
      { NameOfCoef uxn ;
        EntityType NodesOf ; NameOfConstraint Displacement_x ; }
-      If ( FE_Degree == 2 )
+      If ( FE_Order == 2 )
         { NameOfCoef uxn2 ;
 	   EntityType EdgesOf ; NameOfConstraint Displacement_x ; }
     EndIf
@@ -187,7 +187,7 @@ FunctionSpace {
    BasisFunction {
      { Name syn ; NameOfCoef uyn ; Function BF_Node ;
        Support Dom_H_u_Mec ; Entity NodesOf[ All ] ; }
-     If ( FE_Degree == 2 )
+     If ( FE_Order == 2 )
        { Name syn2 ; NameOfCoef uyn2 ; Function BF_Node_2E ;
          Support Dom_H_u_Mec; Entity EdgesOf[ All ] ; }
     EndIf
@@ -195,7 +195,7 @@ FunctionSpace {
    Constraint {
      { NameOfCoef uyn ;
        EntityType NodesOf ; NameOfConstraint Displacement_y ; }
-      If ( FE_Degree == 2 )
+      If ( FE_Order == 2 )
      { NameOfCoef uyn2 ;
        EntityType EdgesOf ; NameOfConstraint Displacement_y ; }
      EndIf
@@ -222,7 +222,7 @@ Jacobian {
 Integration {
  { Name Gauss_v;
    Case {
-      If (FE_Degree == 1)
+      If (FE_Order == 1)
      { Type Gauss;
        Case {
          { GeoElement Line       ; NumberOfPoints  3; }
@@ -319,7 +319,7 @@ PostProcessing {
 PostOperation {
  { Name pos; NameOfPostProcessing Elast_u;
    Operation {
-      If(FE_Degree == 1)
+      If(FE_Order == 1)
 	Print[ sig_xx, OnElementsOf Wrench, File "sigxx.pos" ];
        Print[ u, OnElementsOf Wrench, File "u.pos" ];
      Else

--- a/Electrostatics/microstrip.pro
+++ b/Electrostatics/microstrip.pro
@@ -32,7 +32,7 @@
   We consider here the special case where rho = 0 to model a conducting
   microstrip line on top of a dielectric substrate. A Dirichlet boundary
   condition sets the potential to 1 mV on the boundary of the microstrip line
-   (called "Electrode" below) and 0 V on the ground. A homogeneous Neumann
+   (called "Electrode" below) and to 0 V on the ground. A homogeneous Neumann
   boundary condition (zero flux of the displacement field, i.e. n.d = 0) is
   imposed on the left boundary of the domain to account for the symmetry of the
   problem, as well as on the top and right boundaries that truncate the
@@ -59,8 +59,8 @@ Group {
     Vol_xxx groups contain only volume elements of the mesh (triangles here).
     Sur_xxx groups contain only surface elements of the mesh (lines here).
-     Since there are no non-homogeneous Neumann conditions in the model,
+     Since there are no non-homogeneous Neumann conditions in this particular
-     Sur_Neu_Ele is defined as empty.
+     example, Sur_Neu_Ele is defined as empty.
     */
  Vol_Ele = Region[ {Air, Diel1} ];
@@ -141,18 +141,22 @@ Jacobian {
     the reference elements (defined in standardized unit cells) over which
     integration is performed. "Vol" represents the classical 1-to-1 mapping
     between identical spatial dimensions, i.e. in this case a reference
-     triangle/quadrangle onto triangles/quadrangles in the z=0 plane
+     triangle/quadrangle onto triangles/quadrangles in the z=0 plane (2D <->
-     (2D <-> 2D). "Sur" would be used to map the reference triangle/quadrangle
+     2D). "Sur" is used to map the reference triangle/quadrangle onto
-     onto triangles/quadrangles in a 3D space (2D <-> 3D), or to map the
+     triangles/quadrangles in a 3D space (2D <-> 3D), or to map the reference
-     reference line segment onto segments in 2D space (1D <-> 2D). "Lin" would
+     line segment onto segments in 2D space (1D <-> 2D). "Lin" is used to map
-     be used to map the reference line segment onto segments in 3D space (1D <->
+     the reference line segment onto segments in 3D space (1D <-> 3D). */
-     3D). */
  { Name Vol ;
    Case {
      { Region All ; Jacobian Vol ; }
    }
  }
+  { Name Sur ;
+    Case {
+      { Region All ; Jacobian Sur ; }
+    }
+  }
 }
 Integration {
@@ -160,7 +164,8 @@ Integration {
  { Name Int ;
    Case { { Type Gauss ;
-             Case { { GeoElement Triangle    ; NumberOfPoints  4 ; }
+             Case { { GeoElement Line        ; NumberOfPoints  4 ; }
+                    { GeoElement Triangle    ; NumberOfPoints  4 ; }
                    { GeoElement Quadrangle  ; NumberOfPoints  4 ; } }
      }
    }
@@ -253,6 +258,15 @@ Formulation {
    Equation {
      Integral { [ epsilon[] * Dof{d v} , {d v} ];
 	    In Vol_Ele; Jacobian Vol; Integration Int; }
+   /* Additional Integral terms could be added here. For example, the
+      following term would account for non-homogeneous Neumann boundary
+      conditions, provided that the function nd[] is defined:
+      Integral { [ nd[] , {v} ];
+	    In Sur_Neu_Ele; Jacobian Sur; Integration Int; }
+      All the terms are added, and an implicit "= 0" is considered at the end. */
    }
  }
 }
@@ -294,15 +308,15 @@ PostProcessing {
  { Name EleSta_v; NameOfFormulation Electrostatics_v;
    Quantity {
      { Name v; Value {
-          Local { [ {v} ]; In Dom_Hgrad_v_Ele; Jacobian Vol; }
+          Term { [ {v} ]; In Dom_Hgrad_v_Ele; Jacobian Vol; }
        }
      }
      { Name e; Value {
-          Local { [ -{d v} ]; In Dom_Hgrad_v_Ele; Jacobian Vol; }
+          Term { [ -{d v} ]; In Dom_Hgrad_v_Ele; Jacobian Vol; }
        }
      }
      { Name d; Value {
-          Local { [ -epsilon[] * {d v} ]; In Dom_Hgrad_v_Ele; Jacobian Vol; }
+          Term { [ -epsilon[] * {d v} ]; In Dom_Hgrad_v_Ele; Jacobian Vol; }
        }
      }
    }
@@ -315,8 +329,8 @@ h = 2.e-3; // vertical position of the cut
 PostOperation {
  { Name Map; NameOfPostProcessing EleSta_v;
     Operation {
-       Print [ v, OnElementsOf Dom_Hgrad_v_Ele, File "mStrip_v.pos" ];
+       Print [ v, OnElementsOf Vol_Ele, File "mStrip_v.pos" ];
-       Print [ e, OnElementsOf Dom_Hgrad_v_Ele, File "mStrip_e.pos" ];
+       Print [ e, OnElementsOf Vol_Ele, File "mStrip_e.pos" ];
       Print [ e, OnLine {{e,h,0}{14.e-3,h,0}}{60}, File "Cut_e.pos" ];
     }
  }

--- a/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro
+++ b/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro
@@ -15,7 +15,7 @@ DefineConstant[
  TimeInit = 0, // intial time (for time-domain simulations)
  TimeFinal = 1/50, // final time (for time-domain simulations)
  DeltaTime = 1/500, // time step (for time-domain simulations)
-  InterpolationOrder = 1 // finite element order
+  FE_Order = 1 // finite element order
  Val_Rint = 0, // interior radius of annulus shell transformation region (Vol_Inf_Mag)
  Val_Rext = 0 // exterior radius of annulus shell  transformation region (Vol_Inf_Mag)
 ];
@@ -134,7 +134,7 @@ FunctionSpace {
        Support Dom_Mag; Entity GroupsOfNodesOf[Sur_Perfect_Mag]; }
      // additional basis functions for 2nd order interpolation
-      If(InterpolationOrder == 2)
+      If(FE_Order == 2)
        { Name s_e; NameOfCoef a_e; Function BF_PerpendicularEdge_2E;
          Support Vol_Mag; Entity EdgesOf[All]; }
      EndIf
@@ -150,7 +150,7 @@ FunctionSpace {
      { NameOfCoef I;
        EntityType GroupsOfNodesOf; NameOfConstraint Current_2D; }
-      If(InterpolationOrder == 2)
+      If(FE_Order == 2)
        { NameOfCoef a_e;
          EntityType EdgesOf; NameOfConstraint MagneticVectorPotential_2D_0; }
      EndIf
@@ -386,10 +386,19 @@ PostProcessing {
  { Name MagDyn_a_2D; NameOfFormulation MagDyn_a_2D;
    PostQuantity {
      // In 2D, a is a vector with only a z-component: (0,0,az)
-      { Name a; Value { Term { [ {a} ]; In Vol_Mag; Jacobian Vol; } } }
+      { Name a; Value {
+          Term { [ {a} ]; In Vol_Mag; Jacobian Vol; }
+        }
+      }
      // The equilines of az are field lines (giving the magnetic field direction)
-      { Name az; Value { Term { [ CompZ[{a}] ]; In Vol_Mag; Jacobian Vol; } } }
+      { Name az; Value {
-      { Name b; Value { Term { [ {d a} ]; In Vol_Mag; Jacobian Vol; } } }
+          Term { [ CompZ[{a}] ]; In Vol_Mag; Jacobian Vol; }
+        }
+      }
+      { Name b; Value {
+          Term { [ {d a} ]; In Vol_Mag; Jacobian Vol; }
+        }
+      }
      { Name h; Value {
          Term { [ nu[] * {d a} ]; In Vol_Mag; Jacobian Vol; }
          Term { [ -nu[] * br[] ]; In Vol_M_Mag; Jacobian Vol; }
@@ -401,9 +410,7 @@ PostProcessing {
          Term { [ Vector[0,0,0] ]; In Vol_Mag; Jacobian Vol; } // to force a vector result out of sources
        }
      }
-      { Name j;
+      { Name j; Value {
-        // Only correct in MagDyn
-        Value {
          Term { [ -sigma[] * (Dt[{a}]+{ur}/CoefGeos[]) ]; In Vol_C_Mag; Jacobian Vol; }
          Term { [ js0[] ]; In Vol_S0_Mag; Jacobian Vol; }
          Term { [ (js0[]*Vector[0,0,1])*{ir} ]; In Vol_S_Mag; Jacobian Vol; }
@@ -412,22 +419,17 @@ PostProcessing {
          Term { [ -Ysur[] * (Dt[{a}]+{ur}/CoefGeos[]) ]; In Sur_Imped_Mag; Jacobian Sur; }
        }
      }
+      { Name JouleLosses; Value {
-      { Name JouleLosses;
-        Value {
          Integral { [ CoefPower * sigma[]*SquNorm[Dt[{a}]+{ur}/CoefGeos[]] ];
            In Vol_C_Mag; Jacobian Vol; Integration Gauss_v; }
          Integral { [ CoefPower * 1./sigma[]*SquNorm[js0[]] ];
            In Vol_S0_Mag; Jacobian Vol; Integration Gauss_v; }
 	  Integral { [ CoefPower * 1./sigma[]*SquNorm[(js0[]*Vector[0,0,1])*{ir}] ];
            In Vol_S_Mag; Jacobian Vol; Integration Gauss_v; }
          Integral { [ CoefPower * Ysur[]*SquNorm[Dt[{a}]+{ur}/CoefGeos[]] ];
            In Sur_Imped_Mag; Jacobian Sur; Integration Gauss_v; }
 	}
      }
      { Name U; Value {
          Term { [ {U} ]; In Vol_C_Mag; }
          Term { [ {Us} ]; In Vol_S_Mag; }
@@ -436,7 +438,6 @@ PostProcessing {
          EndIf
        }
      }
      { Name I; Value {
          Term { [ {I} ]; In Vol_C_Mag; }
          Term { [ {Is} ]; In Vol_S_Mag; }
@@ -445,20 +446,28 @@ PostProcessing {
          EndIf
        }
      }
    }
  }
  { Name MagSta_a_2D; NameOfFormulation MagSta_a_2D;
    PostQuantity {
-      // In 2D, a is a vector with only a z-component: (0,0,az)
+      { Name a; Value {
-      { Name a; Value { Term { [ {a} ]; In Vol_Mag; Jacobian Vol; } } }
+          Term { [ {a} ]; In Vol_Mag; Jacobian Vol; }
-      // The equilines of az are field lines (giving the magnetic field direction)
+        }
-      { Name az; Value { Term { [ CompZ[{a}] ]; In Vol_Mag; Jacobian Vol; } } }
+      }
-      { Name b; Value { Term { [ {d a} ]; In Vol_Mag; Jacobian Vol; } } }
+      { Name az; Value {
-      { Name h; Value { Term { [ nu[] * {d a} ]; In Vol_Mag; Jacobian Vol; } } }
+          Term { [ CompZ[{a}] ]; In Vol_Mag; Jacobian Vol; }
-      { Name j;
+        }
-        Value {
+      }
+      { Name b; Value {
+          Term { [ {d a} ]; In Vol_Mag; Jacobian Vol; }
+        }
+      }
+      { Name h; Value {
+          Term { [ nu[] * {d a} ]; In Vol_Mag; Jacobian Vol; }
+        }
+      }
+      { Name j; Value {
          Term { [ js0[] ]; In Vol_S0_Mag; Jacobian Vol; }
          Term { [ (js0[]*Vector[0,0,1])*{ir} ]; In Vol_S_Mag; Jacobian Vol; }
          Term { [ Vector[0,0,0] ]; In Vol_Mag; Jacobian Vol; }

--- a/Magnetostatics/electromagnet.pro
+++ b/Magnetostatics/electromagnet.pro
@@ -42,15 +42,15 @@
   This model computes the static magnetic field produced by a DC current. This
   corresponds to a "magnetostatic" physical model, obtained by combining the
-   time-invariant Maxwell-Ampere equation (Curl h = js, with h the magnetic
+   time-invariant Maxwell-Ampere equation (curl h = js, with h the magnetic
   field and js the source current density) with Gauss' law (Div b = 0, with b
   the magnetic flux density) and the magnetic constitutive law (b = mu h, with
   mu the magnetic permeability).
   Since Div b = 0, b can be derived from a vector magnetic potential a, such
-   that b = Curl a. Plugging this potential in Maxwell-Ampere's law and using
+   that b = curl a. Plugging this potential in Maxwell-Ampere's law and using
   the constitutive law leads to a vector Poisson equation in terms of the
-   magnetic vector potential: Curl(nu Curl a) = js, where nu = 1/mu is
+   magnetic vector potential: curl(nu curl a) = js, where nu = 1/mu is
   the reluctivity. */
 Group {
@@ -74,34 +74,35 @@ Group {
     - Vol_Mag     : full volume domain
     - Vol_S_Mag   : region where the current source js is defined
     - Vol_Inf_Mag : region where the infinite ring geometric transformation is applied
-     - Sur_Dir_Mag : homogeneous Dirichlet part of the model boundary
+     - Sur_Dir_Mag : part of the boundary with homogenous Dirichlet conditions
-     - Sur_Neu_Mag : homogeneous Neumann part of the model boundary
+     - Sur_Neu_Mag : part of the boundary with non-homogeneous Neumann conditions
  */
  Vol_Mag     = Region[ {Air, AirInf, Core, Ind} ];
  Vol_S_Mag   = Region[ Ind ];
  Vol_Inf_Mag = Region[ AirInf ];
  Sur_Dir_Mag = Region[ {Surface_bn0, Surface_Inf} ];
-  Sur_Neu_Mag = Region[ Surface_ht0 ];
+  Sur_Neu_Mag = Region[ {} ]; // empty
 }
 /* The weak formulation for this problem is derived in a similar way to the
   electrostatic weak formulation from Tutorial 1. The main difference is that
-   the fields are now vector-valued, and that we have one linear term in
+   the fields are now vector-valued, and that we have one linear (source) term
-   addition to the bilinear term. The weak formulation reads: find a such that
+   in addition to the bilinear term. The weak formulation reads: find a such
+   that
-   (Curl(nu Curl a), a')_Vol_Mag = (js, a')_Vol_S_Mag
+   (curl(nu curl a), a')_Vol_Mag = (js, a')_Vol_S_Mag
   holds for all test functions a'. After integration by parts it reads: find a
   such that
-   (nu Curl a, Curl a')_Vol_Mag + (n x (nu Curl a), a')_Bnd_Vol_Mag = (js, a')_Vol_S_Mag
+   (nu curl a, curl a')_Vol_Mag + (n x (nu curl a), a')_Bnd_Vol_Mag = (js, a')_Vol_S_Mag
-   In this electromagnet model the second (boundary term) vanishes, as there is
+   In this electromagnet model the second (boundary) term vanishes, as there is
-   either no test function a' (on the Dirichlet boundary), or "n x (nu Curl a) =
+   either no test function a' (on the Dirichlet boundary), or "n x (nu curl a) =
   n x h" is zero (on the homogeneous Neumann boundary). We are thus eventually
   looking for functions a such that
-   (nu Curl a, Curl a')_Vol_Mag = (js, a')_Vol_S_Mag
+   (nu curl a, curl a')_Vol_Mag = (js, a')_Vol_S_Mag
   holds for all a'. */
@@ -228,6 +229,8 @@ Formulation {
      { Name js; Type Local; NameOfSpace Hregion_j_Mag_2D; }
    }
    Equation {
+      // all terms on the left-hand side (hence the "-" sign in front of
+      // Dof{js}):
      Integral { [ nu[] * Dof{d a} , {d a} ];
        In Vol_Mag; Jacobian Vol; Integration Int; }
      Integral { [ -Dof{js} , {a} ];
@@ -252,27 +255,27 @@ PostProcessing {
    Quantity {
      { Name a;
        Value {
-          Local { [ {a} ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
+          Term { [ {a} ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
        }
      }
      { Name az;
        Value {
-          Local { [ CompZ[{a}] ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
+          Term { [ CompZ[{a}] ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
        }
      }
      { Name b;
        Value {
-          Local { [ {d a} ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
+          Term { [ {d a} ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
        }
      }
      { Name h;
        Value {
-          Local { [ nu[] * {d a} ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
+          Term { [ nu[] * {d a} ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
        }
      }
      { Name js;
        Value {
-          Local { [ {js} ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
+          Term { [ {js} ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
        }
      }
    }

--- a/PotentialFlow/magnus.pro
+++ b/PotentialFlow/magnus.pro
@@ -315,29 +315,30 @@ PostProcessing{
          Term{ [ argVTrail[{d phi}] ]; In Dom_Vh; Jacobian Vol; }
        }
      }
+      { Name Circ; Value {
-      { Name Circ; Value { Term { [ {Circ} ]; In Sur_Cut; } } }
+          Term { [ {Circ} ]; In Sur_Cut; }
+        }
-      { Name Dmdt; Value { Term { [ {Dmdt} ]; In Sur_Cut; } } }
+      }
+      { Name Dmdt; Value {
+          Term { [ {Dmdt} ]; In Sur_Cut; }
+        }
+      }
      // Kutta-Jukowski approximation for Lift
-      { Name LiftKJ; Value { Term { [ -rho[]*{Circ}*Velocity ]; In Sur_Cut; } } }
+      { Name LiftKJ; Value {
+          Term { [ -rho[]*{Circ}*Velocity ]; In Sur_Cut; }
+        }
+      }
      // Lift computed with the real pressure field
-      { Name Lift;
+      { Name Lift; Value {
-         Value {
          Integral { [ -0.5*rho[]*SquNorm[{d phi}]*CompY[Normal[]] ];
            In Dom_Vh ; Jacobian Sur ; Integration Int; }
        }
      }
+      { Name circulation; Value {
-      { Name circulation;
-         Value {
          Integral { [ {d phi} * Tangent[] ] ;
            In Dom_Vh ; Jacobian Sur  ; Integration Int; }
        }
      }
    }
  }
 }