[Getdp] Concatenation of field operators

[Thomas_W]

Hi Michael,

we had u as discretized quantity, and wanted to evaluate |grad|grad u||.

I have tried your suggestion with the new quantity. I attach the  
files. The result is that the computed result is just zero everywhere.

Are there some other ideas how to compute the term?  - Probably, there  
is still a mistake in the files.

Thanks, in any event!
Thomas


Zitat von michael.asam at infineon.com:

> Hi Thomas,
>
> I think this cannot be done in the post-processing.
> Probably one solution would be to calculate the first gradient
> in the Formulation{} part by adding a new quantity (e.g. dv)
> and do a kind of protection:
>
> Formulation{
> ...
>   Galerkin { [ -Norm[{v}], {dv} ]; In All; Integration I1; Jacobian JVol;  }
>   Galerkin { [ Dof{dv}, {dv} ]; In All; Integration I1; Jacobian JVol;  }
> }
>
> You have to make sure that v has at least 2nd order basis functions
> in order to be twice differentiable.
>
> I hope this works :-)
>
> Best regards
> Michael
>
>
>
>
>
> -----Original Message-----
> From: getdp-bounces at ace20.montefiore.ulg.ac.be  
> [mailto:getdp-bounces at ace20.montefiore.ulg.ac.be] On Behalf Of  
> Thomas_W
> Sent: Monday, October 15, 2012 3:03 PM
> To: getdp at geuz.org
> Subject: [Getdp] Concatenation of field operators
>
> Dear all,
>
> I try to evaluate the following quantity in the Post-processing part,
> where v is the discretized quantity:
> |Grad (|Grad (v)|)|.
>
> I tried variants such as:
>
> PostProcessing {
>    { Name The; NameOfFormulation f_v;
>      Quantity {
>        { Name jj;
>          Value{
>          Local{ [ {d  {d v}*{d v} } ] ; In All;Jacobian JVol; }
>          }
>        }
>      }
>    }
> }
>
> or also
> Local{ [ {d  {Norm[{v}]} } ] ; In All;Jacobian JVol; }.
>
> All these variants yield an error, namely 'syntax error ({)' in the
> corresponding line. It probably has to do with the fact that only v,
> and not Norm[{v}] is the original discretized quantity.
>
> Do you have a solution how to evaluate such a function with GetDP?
>
> Thank you,
> Thomas
>
>
> _______________________________________________
> getdp mailing list
> getdp at geuz.org
> http://www.geuz.org/mailman/listinfo/getdp
>
>


-------------- next part --------------
// File "LaplacianNeumann.geo".

// We include the file containing the numbering of the geometry.
// This is usefull at the end of this file, and used to "synchronise" GMSH and GetDP
Include "param.geo";

//Caracteristic length of the finite elements (reffinement is also possible after the mesh is built):
lc = 0.05;
// This parameter could be placed for instance in "param.geo", to separate more easyly the geometry
// 	and the discretization parameters.

// The parameters of the border of the domain :
x_max = 1;
x_min = 0;
y_max= 1;
y_min = 0;

//Creation of the 4 angle points of the domain Omega (=square)
p1 = newp; Point(p1) = {x_min,y_min,0,lc};
p2 = newp; Point(p2) = {x_min,y_max,0,lc};
p3 = newp; Point(p3) = {x_max,y_max,0,lc};
p4 = newp; Point(p4) = {x_max,y_min,0,lc};
// Remarks:
// -"newp" is a GMSH function that give the first available number for describing a point.
// 	For any other entity, like Line, Surface, etc. We recommand the use of "newreg" (see below).
// - By default, GMSH create a 3D domain. The z-coordinate must always be precised.


//The four edges of the square
L1 = newreg; Line(L1) = {p1,p2};
L2 = newreg; Line(L2) = {p2,p3};
L3 = newreg; Line(L3) = {p3,p4};
L4 = newreg; Line(L4) = {p4,p1};

// Line Loop (= boundary of the square)
Bound = newreg; Line Loop(Bound) = {L1,L2,L3,L4};

//Surface of the square
SurfaceOmega = newreg; Plane Surface(SurfaceOmega) = {Bound};

// To conclude, we define the physical entities, that is "what GetDP could see/use".
// "Omega" is a number imported from the file "param.geo".
Physical Surface(Omega) = {SurfaceOmega};
// Do not forget to let a blank line at the end, this could make GMSH crash...
Physical Line(Left)={L1}; 
Physical Line(Top) ={L2};
Physical Line(Right)={L3};
Physical Line(Bottom)={L4};


-------------- next part --------------

Include "param.geo";

Group{
Omega = Region[{Omega}];
Top =#Top;
Bottom = #Bottom;
Left = #Left;
Right = #Right;
Boundlined = Region[{Left,Right,Top,Bottom}];
}

Function{
Coef = Pi;
f[] = (1+2*Coef*Coef)*Cos[Coef*X[]]*Cos[Coef*Y[]];
}



Jacobian {
  { Name JVol ;
    Case {
      { Region All ; Jacobian Vol ; }
    }
  }
  { Name JSur ;
    Case {
      { Region All ; Jacobian Sur ; }
    }
  }
  { Name JLin ;
    Case {
      { Region All ; Jacobian Lin ; }
    }
  }
}


Integration {
  { Name I1 ;
    Case {
      { Type Gauss ;
        Case {
          { GeoElement Point ; NumberOfPoints  1 ; }
          { GeoElement Line ; NumberOfPoints  4 ; }
          { GeoElement Triangle ; NumberOfPoints  6 ; }
          { GeoElement Quadrangle ; NumberOfPoints 7 ; }
          { GeoElement Tetrahedron ; NumberOfPoints 15 ; }
          { GeoElement Hexahedron ; NumberOfPoints 34 ; }
        }
      }
    }
  }
}




Constraint {
  { Name Sta_V ; Type Assign;
    Case {
      { Region Bottom; Value 0; }
      { Region Top; Value 1;}
    }
    Name Sta_V ; Type Assign;
    Case {
      { Region Boundlined; Value 0; }
    }
  }
}


FunctionSpace{ // 2nd order basis function 
  { Name Vh; Type Form0;
    BasisFunction{
    {Name wn; NameOfCoef vn; Function BF_Node;    Support Omega; Entity NodesOf[All];}
    {Name wn2; NameOfCoef vn2; Function BF_Node_2E;    Support Omega; Entity EdgesOf[All];}
    }
    
    Constraint {
     { NameOfCoef vn2; EntityType EdgesOf ; NameOfConstraint Sta_V; }
    }
  }

    { Name Vh2; Type Form0;
    BasisFunction{
      {Name wn; NameOfCoef vnr; Function BF_Node;    Support Omega; Entity NodesOf[All];}
    } 
  }
}



Formulation{
  {Name Laplacian; Type FemEquation;
    Quantity{
      {Name u; Type Local; NameOfSpace Vh;}
      {Name du; Type Local; NameOfSpace Vh2;} // this is for the evaluation of the change of current density
    }
   Equation{
      Galerkin{ [Dof{Grad u}, {Grad u}];
        In Omega; Jacobian JVol; Integration I1;}

      Galerkin{ [  Dof{u}, {u}];
       In Omega; Jacobian JVol; Integration I1;}
     
      Galerkin{ [-f[], {u}];
        In Omega; Jacobian JVol; Integration I1;}
        
        
        
      Galerkin{ [  Dof{du}, {du}];
       In Omega; Jacobian JVol; Integration I1;}
       
       
      Galerkin{ [  -Norm[{Grad u}], {du}];
       In Omega; Jacobian JVol; Integration I1;}
        
        
    }
  }
}

Resolution{
  {Name Laplacian;
    System{
      {Name Syst; NameOfFormulation Laplacian;}
    }
    Operation{
      Generate[Syst]; Solve[Syst]; SaveSolution[Syst];
    }
  }
}


PostProcessing{
  {Name Laplacian; NameOfFormulation Laplacian;
    Quantity{
      {Name u; Value {Local{[{u}];In Omega;Jacobian JVol;}}}
      // this gives the current density vector grad u
      {Name CDe1; Value {Local{[{Grad u}]; In Omega; Jacobian JVol;}}}   
      // this does not give the absolute value |grad u| of the current density vector, but zero   
      {Name CDe2 ; Value { Term { [{du}]   ; In Omega ; Jacobian JVol;} } } 
      // this does not give |grad|grad u||
      {Name CDe3 ; Value { Term { [Norm[{d du}]]   ; In Omega ; Jacobian JVol;} } } 
    
    }
  }
}

PostOperation{
  {Name Map_u; NameOfPostProcessing Laplacian;
    Operation{
      Print[u, OnElementsOf Omega, File "u.pos"];
      Print[CDe1, OnElementsOf Omega, File "CDe1.pos"];
      Print[CDe2, OnElementsOf Omega, File "CDe2.pos"];
      Print[CDe3, OnElementsOf Omega, File "CDe3.pos"];
    }
  }
}

-------------- next part --------------
// File "param.geo"

//Numbers that caracterise the interior of the square (Omega) and its boundary (Gama):
Omega = 1000;
Top = 111;
Left = 112;
Bottom = 113;
Right = 114;

// Three remarks on these numbers :
// - They are arbitrary choosen.
// - They are placed in a separated file to be readable by both GMSH and GetDP.
// - "Gamma" is a special word used by GMSH/GetDP, that is why the boundary is named "Gama", with one "m"...
// Do not forget to let a blank line at the end, this could make GMSH crash...