Determining the direction of the normal between two regions.

On 20 July 2011 15:56, Evan Lezar <email address hidden> wrote:
> New question #165463 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/165463
>
> Hi,
>
> Consider the following:
>
> A mesh consisting of two regions/domains, A and B, with A contained in B (A is a SubDomain of B). How would one go about determining the orientation of the normal vector pointing from A to B when integrating over the exterior facets of A (the surface separating the two regions)?

As far as I know, there's no easy way of doing this.

> These integrals are required for the computation of flux quantities and the near to far-field transform (NTFF) in the electromagnetic problems that we are trying to solve.
>
> For the example given below, consisting of two concentric square domains, the integration results in
> left: 0.5
> top: 0.0
> right: 0.5
> bottom: 0.0,
> where -0.5, 0.0, 0.5, 0.0 is the desired result -- indicating that the normal points outward from the inner square domain.
>
> Example:
> """construct a 1x1 mesh that has inside it a 0.5x0.5 square mesh.
>
> Integrate the dot product of the facet normal with a unit field in the x direction along each one of the
> sides of the interior square (0.5x0.5). The integral of the two vertical sides should equal the length of each side,
> with the sign indicating the direction of the normal. For the normal pointing out of the inner square,
> the left (100) integral should be -0.5 and the right (102) integral should be +0.5.
> """
> from dolfin import *
>
> class Left ( SubDomain ):
>    def inside (self, x, on_boundary ):
>        if x[0] == 0.25 and ( x[1] >= 0.25 and x[1] <= 0.75 ):
>            return True
>        return False
>
> class Right ( SubDomain ):
>    def inside (self, x, on_boundary ):
>        if x[0] == 0.75 and ( x[1] >= 0.25 and x[1] <= 0.75 ):
>            return True
>        return False
>
> class Top ( SubDomain ):
>    def inside (self, x, on_boundary ):
>        if x[1] == 0.75 and ( x[0] >= 0.25 and x[0] <= 0.75 ):
>            return True
>        return False
>
> class Bottom ( SubDomain ):
>    def inside (self, x, on_boundary ):
>        if x[1] == 0.25 and ( x[0] >= 0.25 and x[0] <= 0.75 ):
>            return True
>        return False
>
> sides = ['left', 'top', 'right', 'bottom']
>
> mesh = UnitSquare ( 8, 8 )
> V = FunctionSpace(mesh, "CG", 1)
> u = TrialFunction(V)
> v = TestFunction(V)
> # a uniform x-directed field
> E = Expression ( ("1.0", "0.0") )
>
> interior_edges = MeshFunction ( 'uint', mesh, 1 )
> interior_edges.set_all ( 0 )
>
> left_side = Left ( )
> left_side.mark( interior_edges, 100 )
> top_side = Top ( )
> top_side.mark( interior_edges, 101 )
> right_side = Right ( )
> right_side.mark( interior_edges, 102 )
> bottom_side = Bottom ( )
> bottom_side.mark( interior_edges, 103 )
>
> for interior_edge in [100, 101, 102, 103]:
>    n = V.cell().n
>    l_p = dot( n('+'), E('+'))*dS(interior_edge)
>    print sides[interior_edge-100], assemble ( l_p, interior_facet_domains=interior_edges, mesh=mesh )
>
> """End code snippet"""

Not too familiar with the Python interface, but you should be able to
do something like:

C = Constant("triangle")

Fill C with 1.0 for triangles inside your subdomain, and 0.0 outside
(don't know the correct syntax), and then
define your l_p:

l_p = (C('+')dot( n('+'), E('+') + C('-')dot( n('-'), E('-'))*dS(interior_edge)

that should do the trick.

Kristian

> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>