From 7c8956f7ae55c1a04bf25d3ccf2b9d5f65bca410 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Fri, 16 Mar 2018 12:05:05 +0100
Subject: [PATCH] up
---
Elasticity/wrench2D.pro | 14 ++---
Electrostatics/microstrip.pro | 44 ++++++++++-----
Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro | 61 ++++++++++++---------
Magnetostatics/electromagnet.pro | 39 +++++++------
PotentialFlow/magnus.pro | 39 ++++++-------
5 files changed, 112 insertions(+), 85 deletions(-)
diff --git a/Elasticity/wrench2D.pro b/Elasticity/wrench2D.pro
index a1f8dc1..d7a27a0 100644
--- 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
diff --git a/Electrostatics/microstrip.pro b/Electrostatics/microstrip.pro
index 43bf1da..3e48414 100644
--- 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,
- Sur_Neu_Ele is defined as empty.
+ Since there are no non-homogeneous Neumann conditions in this particular
+ 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
- (2D <-> 2D). "Sur" would be used to map the reference triangle/quadrangle
- onto triangles/quadrangles in a 3D space (2D <-> 3D), or to map the
- reference line segment onto segments in 2D space (1D <-> 2D). "Lin" would
- be used to map the reference line segment onto segments in 3D space (1D <->
- 3D). */
+ triangle/quadrangle onto triangles/quadrangles in the z=0 plane (2D <->
+ 2D). "Sur" is used to map the reference triangle/quadrangle onto
+ triangles/quadrangles in a 3D space (2D <-> 3D), or to map the reference
+ line segment onto segments in 2D space (1D <-> 2D). "Lin" is used to map
+ the reference line segment onto segments in 3D space (1D <-> 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 [ e, OnElementsOf Dom_Hgrad_v_Ele, File "mStrip_e.pos" ];
+ Print [ v, OnElementsOf Vol_Ele, File "mStrip_v.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" ];
}
}
diff --git a/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro b/Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro
index 145c95a..0ec6524 100644
--- 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 b; Value { Term { [ {d a} ]; In Vol_Mag; Jacobian Vol; } } }
+ { Name az; Value {
+ 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,33 +410,26 @@ PostProcessing {
Term { [ Vector[0,0,0] ]; In Vol_Mag; Jacobian Vol; } // to force a vector result out of sources
}
}
- { Name j;
- // Only correct in MagDyn
- Value {
+ { Name j; 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; }
+ Term { [ (js0[]*Vector[0,0,1])*{ir} ]; In Vol_S_Mag; Jacobian Vol; }
Term { [ Vector[0,0,0] ]; In Vol_Mag; Jacobian Vol; }
// Current density in A/m
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 { 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 h; Value { Term { [ nu[] * {d a} ]; In Vol_Mag; Jacobian Vol; } } }
- { Name j;
- Value {
+ { Name a; Value {
+ Term { [ {a} ]; In Vol_Mag; Jacobian Vol; }
+ }
+ }
+ { Name az; Value {
+ 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; }
+ }
+ }
+ { 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; }
diff --git a/Magnetostatics/electromagnet.pro b/Magnetostatics/electromagnet.pro
index d1956f3..0e4126c 100644
--- 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_Neu_Mag : homogeneous Neumann part of the model boundary
+ - Sur_Dir_Mag : part of the boundary with homogenous Dirichlet conditions
+ - 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
- addition to the bilinear term. The weak formulation reads: find a such that
+ the fields are now vector-valued, and that we have one linear (source) term
+ 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
- either no test function a' (on the Dirichlet boundary), or "n x (nu Curl a) =
+ 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) =
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; }
}
}
}
diff --git a/PotentialFlow/magnus.pro b/PotentialFlow/magnus.pro
index eb0252e..77535d0 100644
--- a/PotentialFlow/magnus.pro
+++ b/PotentialFlow/magnus.pro
@@ -226,7 +226,7 @@ Formulation{
Resolution{
{Name PotentialFlow;
System{
- {Name A; NameOfFormulation PotentialFlow;}
+ { Name A; NameOfFormulation PotentialFlow; }
}
Operation{
InitSolution[A];
@@ -287,7 +287,7 @@ Resolution{
PostProcessing{
{ Name PotentialFlow; NameOfFormulation PotentialFlow;
Quantity{
- {Name phi; Value {
+ { Name phi; Value {
Term{ [ {phi} ] ; In Dom_Vh; Jacobian Vol; }
}
}
@@ -299,45 +299,46 @@ PostProcessing{
Term { [ { phiDisc } ] ; In Dom_Vh ; Jacobian Vol ; }
}
}
- {Name velocity; Value {
+ { Name velocity; Value {
Term { [ {d phi} ]; In Dom_Vh; Jacobian Vol; }
}
}
- {Name normVelocity; Value {
+ { Name normVelocity; Value {
Term { [ Norm[{d phi}] ]; In Dom_Vh; Jacobian Vol; }
}
}
- {Name pressure; Value {
+ { Name pressure; Value {
Term { [-0.5*rho[]*SquNorm[ {d phi} ]]; In Dom_Vh; Jacobian Vol; }
}
}
- {Name Angle; Value {
+ { Name Angle; Value {
Term{ [ argVTrail[{d phi}] ]; In Dom_Vh; Jacobian Vol; }
}
}
-
- { Name Circ; Value { Term { [ {Circ} ]; In Sur_Cut; } } }
-
- { Name Dmdt; Value { Term { [ {Dmdt} ]; In Sur_Cut; } } }
-
+ { Name Circ; Value {
+ Term { [ {Circ} ]; 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;
- Value {
+ { Name Lift; 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; }
}
}
-
}
}
}
--
GitLab