[Getdp] How to calculate a numerical quantity with GetDP ?
Christophe Geuzaine
cag32 at case.edu
Tue Sep 19 04:48:57 CEST 2006
Olivier Castany wrote:
> Hello,
>
> after I have solved a problem with GetDP, I would like to compute some
> numerical quantities with GetDP
>
> 1) is it possible ?
Yes, although the computations of some traces can be tricky due to the
way mixed finite elements interpolate the discrete quantities (see next
question).
>
> For example, if I've got a field {T}, how can I do the following :
>
> 2) calculate the flux "- lambda[] * {d T}" throught a surface S of
> the mesh
If {T} is interpolated using nodal elements, {d T} on any surface will
be the tangential gradient. If you want the normal gradient, which is
discontinuous across the surface, you could use the "Trace" operator
(which requires you to specify on which side of the surface you want to
evaluate it), but it's still experimental and not really documented...
Most of the time when need such quantities it's better to define a
"global" quantity (see next question). If you don't want to do that, I
would advise for now to just export the volume data into some
script/program and do the evaluation outside of getdp.
>
> 3) calculate the same flux, but through a surface of the _dual_ mesh,
> close to a surface S of the mesh (how can one define such a
> surface in GetDP ?)
In getdp we compute these global quantities by defining them in a volume
layer on one side of the surface. It is equivalent to using the surface
of the dual mesh. See e.g. my thesis for the details.
>
> 4) calculate the average temperature on a surface S
Try the SurfaceArea[] function.
>
> Following is a sample file for which I would like to compute the
> above-mentioned quantities. The surfaces of interest could be
> surfaces SA and SV.
>
> And a last question (not related to the former) :
>
> 5) In the problem I solved, I need to write "Support Region [{D,SA}]" in
> the FunctionSpace. Simply writing "Support D" yields irrelevant results.
> However, all the nodes of SA are nodes of D.
> Could somebody explain why I need to do that ? What exaclty is this
> "Support group-def" ? Isn't it the set of nodes in the "group-def" ?
>
Yes: even if the nodes are contained in both groups, the support of the
basis functions is what is considered when we build the equations. So
you need to specify DA explicitly in the support.
> Have a nice day,
>
> Olivier C.
>
> -------------------
> barre.geo :
>
> /*
> Metal rod
> */
>
> lc_ext = 0.03;
>
> rayon = 0.05;
> longueur = 0.5;
>
> Point(1) = {0,0,0,1};// longueur caractéristique du centre sans utilité
> Point(2) = {rayon,0,0,lc_ext};
> Point(3) = {-rayon,0,0,lc_ext};
> Point(4) = {0,-rayon,0,lc_ext};
> Point(5) = {0,rayon,0,lc_ext};
>
> Circle(1) = {2,1,5};
> Circle(2) = {5,1,3};
> Circle(3) = {3,1,4};
> Circle(4) = {4,1,2};
>
> Line Loop(5) = {1,2,3,4};
> Plane Surface(6) = {5};
>
> Extrude {0,0,longueur} {
> Surface{6}; Layers{ 30,1 }; Recombine;
> }
>
> Physical Surface(31) = {6};
> Physical Surface(32) = {27,15,19,23,28};
> Physical Volume(33) = {1};
>
> /* Remark : gmsh statistics show :
>
> Geometry : 24 lines... it should be : 12
> Mesh : 296 nodes on lines... it should be : 132
> 400 nodes on surfaces... it should be : 260
> */
>
> -------------
> barre.pro :
>
> /*
> Temperature in a metal rod :
>
> - the temperature is imposed at one end (surface SV)
> - convection happens on all other surfaces (surface SA) (j.n = h(T -
> T_ext))
> - conduction of heat in the material (volume D) (j = - lambda Grad(T))
> - equilibrium situation (Div(j)=0)
> */
>
> Group{
> D = Region[33];// Entire volume
> SV = Region[31];// Surface with imposed temperature
> SA = Region[32];// Other surfaces with air convection
> }
>
>
> Function {
> h[] = 30.;// Wm-2K-1 (strong natural air convection)
> lambda[] = 20.;// Wm-1K-1 (some highly allied steel)
> T_ext[] = 30.;// °C (air temperature)
> }
>
>
> // Temperature imposed on surface SV :200°C
> Constraint {
> { Name temperature_sur_SV ; Type Assign;
> Case {
> { Region SV ; Value 200. ; }
> }
> }
> }
>
> FunctionSpace {
> { Name temperature_dans_barre ; Type Form0 ;
> BasisFunction {
> { Name w_n ; NameOfCoef T_n ; Function BF_Node ;
> Support Region [{D,SA}] ; Entity NodesOf[ All ] ; }
> }
> Constraint {
> { NameOfCoef T_n ; EntityType NodesOf ;
> NameOfConstraint temperature_sur_SV ; }
> }
> }
> }
>
> Jacobian {
> { Name JVol ;
> Case {
> { Region All ; Jacobian Vol ; }
> }
> }
> { Name JSur ;
> Case {
> { Region All ; Jacobian Sur ; }
> }
> }
> }
>
> Integration {
> { Name I ;
> Case {
> { 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 ; }
> }
> }
> }
> }
> }
>
> /****************** Integral formulation to solve :
> find T such that :
>
> \int_D lambda Grad(T).Grad(T') + \int_SA h (T - T_ext) T' = 0
> holds for all T' with T' = 0 on SV
> ******************************************/
>
> Formulation {
> { Name formulation_du_probleme ; Type FemEquation ;
> Quantity {
> { Name T ; Type Local ; NameOfSpace temperature_dans_barre ; }
> }
> Equation {
> Galerkin { [ lambda[] * Dof{d T} , {d T} ] ;
> In D ; Jacobian JVol ; Integration I ; }
>
> Galerkin { [ h[] * Dof{T} , {T} ] ;
> In SA ; Jacobian JSur ; Integration I ; }
>
> Galerkin { [ -h[] * T_ext[] , {T} ] ;
> In SA ; Jacobian JSur ; Integration I ; }
>
> }
> }
> }
>
> Resolution {
> { Name resolution_du_probleme ;
> System {
> { Name A ; NameOfFormulation formulation_du_probleme ; }
> }
> Operation {
> Generate[A] ; Solve[A] ; SaveSolution[A] ;
> }
> }
> }
>
> PostProcessing {
> { Name post_processing ; NameOfFormulation formulation_du_probleme ;
> Quantity {
> { Name T; Value { Local { [ {T} ] ; In SA ; Jacobian JSur ; }}}
> }
> }
> }
>
> PostOperation {
> { Name post_operations ; NameOfPostProcessing post_processing;
> Operation {
> Print[ T, OnElementsOf SA, File "T.pos"] ;
> }
> }
> }
>
> ------------------
>
>
> _______________________________________________
> getdp mailing list
> getdp at geuz.org
> http://www.geuz.org/mailman/listinfo/getdp
>
--
Christophe Geuzaine
Assistant Professor, Case Western Reserve University, Mathematics
http://www.case.edu/artsci/math/geuzaine