diff --git a/Elasticity/wrench2D.pro b/Elasticity/wrench2D.pro index a1f8dc128ec65847db1bbf821b4e924a4d4fb221..d7a27a0a68315a1c1fc8881765041cb3b76371b0 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 43bf1dad6e653db345aa0bf64607ccfaf5d73fd3..3e484141a157c359808316083181d14e97b8357a 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 145c95a2facbb957f2ed9655637d15d693c09e8c..0ec6524e312b8ffa4bf4f5c3ef220bbc7f0bc1bb 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 d1956f3fec34051e1d6199ad1e31357bfd8ba445..0e4126c15e2d85e28a5a5ab5564fe863c0d79eb2 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 eb0252e54629079f30b7923c3eb8f00d6ff9b438..77535d0a5fdbaea381cbceebfa2520367a249fc7 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; } } } - } } }