[Getdp] Fluxes through partial boundaries

[Tobias Nähring]

Hello,
I've tried to calculate the fluxes through two parts of a boundary with 
the help of global quantities.

The simple sample problem:
Laplacian equation in (0,1)^2.

Neumann's (natural) boundary conditions at x=0 and x=1,
Dirichlet's boundary conditions at y=0 and y=1:

v(x,0) = 1;
v(x,1) = 0;

I tried to compute the fluxes through the surfaces y=0, 0 < x < 0.5 and 
y=0, 0.5 < x < 1 with the problem definition below.
The problem is that I always get a peak of the potential at the 
connection point (at x=0.5) of the two parts of the lower boundary line.
Even, if I explicitely subtract one half by the other.

Do I choose the wrong ansatz functions?  From my considerations, it 
should work if the two group of nodes for the partial boundaries are 
disjunkt, shouldn't it?

With best regards,
Tobias Naehring

//////////////////////////////////////////////////////////////////////
// Geometry file:
a1=0.5; //< Length of division of surfaces for flux computation
a = 1; //< Side length in x-Direction
b = 1.0; //< Side length in y-direction

fine = 0.1;

Point(1) = {0,0,0,fine};
Point(2) = {a1,0,0,fine};
Point(3) = {a,0,0,fine};
Point(4) = {a,b,0,fine};
Point(5) = {0,b,0,fine};

For i In {1:4}
  Line(i) = {i,i+1};
EndFor

Line(5) = {5,1};

Line Loop (7) = {1:5};

Plane Surface (8) = {7};

Physical Surface (1) = {8};

Physical Line(10) = {1}; ///< 1st partial boundary where we want to know 
the flux
Physical Line(11) = {2}; ///< 2nd partial boundary where we want to know 
the flux
Physical Line(12) = {4}; ///< Plate boundary.

//////////////////////////////////////////////////////////////////////
//Problem definition file:

Group {
    dom = Region[1];
    bdGnd1 = Region[10];
    bdGnd2 = Region[11];
    bdGnd = Region[{bdGnd1,bdGnd2}];
    bd = Region[12];
}

Constraint {
    {
    Name Bc; Type Assign;
    Case
        {
        { Region bd; Value 0.0; }
        }
    }
    {
    Name BcGnd; Type Assign;
    Case
        {
        { Region bdGnd; Value 1.0; } // We need potential 1 on ground 
for the flux computation.
        }
    }
}

FunctionSpace {
    {
    Name funSp;  Type Form0;
    BasisFunction {
        {
        Name bf;
        NameOfCoef v;
        Function BF_Node;
        Support dom;
        Entity NodesOf[All, Not bdGnd];
        }
        {
        Name bfGnd1;
        NameOfCoef vGnd1;
        Function BF_GroupOfNodes;
        Support dom;
        Entity GroupsOfNodesOf[bdGnd1]; //< Here I also tried [bdGnd1, 
Not bdGnd2] with no luck.
        }
        {   
        Name bfGnd2;
        NameOfCoef vGnd2;
        Function BF_GroupOfNodes;
        Support dom;
        Entity GroupsOfNodesOf[bdGnd2]; //< Here I also tried [bdGnd2, 
Not bdGnd1] with no luck.
        }
    }
    GlobalQuantity {
        {
        Name vGlobGnd1;
        Type AliasOf;
        NameOfCoef vGnd1;
        }
        {
        Name fluxGlobGnd1;
        Type AssociatedWith;
        NameOfCoef vGnd1;
        }
        {
        Name vGlobGnd2;
        Type AliasOf;
        NameOfCoef vGnd2;
        }
        {
        Name fluxGlobGnd2;
        Type AssociatedWith;
        NameOfCoef vGnd2;
        }
    }
    Constraint {
        {
        NameOfCoef v;
        EntityType NodesOf;
        NameOfConstraint Bc;
        }
        {
        NameOfCoef vGnd1;
        EntityType GroupsOfNodesOf;
        NameOfConstraint BcGnd;
        }
        {
        NameOfCoef vGnd2;
        EntityType GroupsOfNodesOf;
        NameOfConstraint BcGnd;
        }
    }
    }
}

// Standard choice from the documentation:
Integration {
    {
    Name Int_1;
    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 ; }
            }
        }
    }
    }
}

Jacobian
{
    {
    Name Jac;
    Case {{Region All; Jacobian Sur;}}
    }
}

Formulation {
    {
    Name form; Type FemEquation;
    Quantity {
        {
        Name v; Type Local;
        NameOfSpace funSp;
        }
        {
        Name vDiscreteGnd1; Type Global;
        NameOfSpace funSp[vGlobGnd1];
        }
        {
        Name fluxDiscreteGnd1; Type Global;
        NameOfSpace funSp[fluxGlobGnd1];
        }
        {
        Name vDiscreteGnd2; Type Global;
        NameOfSpace funSp[vGlobGnd2];
        }
        {
        Name fluxDiscreteGnd2; Type Global;
        NameOfSpace funSp[fluxGlobGnd2];
        }
    }
    Equation
        {
        Galerkin
            {
            [ Dof{Grad v} , {Grad v} ];
            In dom;
            Jacobian Jac;
            Integration Int_1;
            }
        GlobalTerm
            {
            [-Dof{fluxDiscreteGnd1}, {vDiscreteGnd1}];
            In bdGnd1;
            }
        GlobalTerm
            {
            [-Dof{fluxDiscreteGnd2}, {vDiscreteGnd2}];
            In bdGnd2;
            }
        }
    }
}

Resolution
{
    {
    Name runSolver;
    System
        {
        {Name sys; NameOfFormulation form;}
        }
    Operation
        { Generate[sys]; Solve[sys]; SaveSolution[sys];}
    }
}

PostProcessing
{
    // The first one is just for verification of the homogeneous field:
    {
    Name postProPotField; NameOfFormulation form;
    Quantity
        {
        {Name v; Value {Local { [{v}]; In dom; Jacobian Jac; }}}
        }
    }
    // THIS IS THE INTERESTING ONE:
    {
    Name postProFlux; NameOfFormulation form;
    Quantity
        {
        {Name flux1; Value {Term {[{fluxDiscreteGnd1}]; In bdGnd1; }}}
        {Name flux2; Value {Term {[{fluxDiscreteGnd2}]; In bdGnd2; }}}
        }
    }
}

PostOperation {
    // The first one is just for checking:
    {
    Name postOpGmsh;
    NameOfPostProcessing postProPotField;
    Operation
        {
        Print[ v, OnElementsOf dom, File "homogen.pos" ];
        }
    }
    // THIS IS THE INTERESTING ONE:
    {
    Name postOpFlux;
    NameOfPostProcessing postProFlux;
    Operation
        {
        Print[ flux1, OnRegion bdGnd1, Format Table, File 
"homogen-flux1.pos" ];
        Print[ flux2, OnRegion bdGnd2, Format Table, File 
"homogen-flux2.pos" ];
        }
    }
}
// End of file homogen.pro
//////////////////////////////////////////////////////////////////////