FEniCS Project

Sum of integrals over each triangle

Asked by Claudia Totzeck

Hi,

I'm trying to solve the continuity equations of the Drift Diffusion model written in the Brezzi, Marini, Pietra Formulation with FEniCS. In the variational formulation occur the sum of integrals over each triangle of the mesh AND the sum of integrals over the boundary of each triangle of the mesh in normal direction. Can anybode give me a hint how to realize this ?

I was thinking about defining each triangle as subdomain... but I believe this is a really messy solution.

Thanks for your help in advance!

Claudia

Question information

Language:
English Edit question
Status:
Solved
For:
FEniCS Project Edit question
Assignee:
No assignee Edit question
Solved by:
Marie Rognes
Solved:
Last query:
Last reply:
Revision history for this message
David Ham (david-ham) said : 		#1

I'm not familiar with that equation set, but sums of integrals over cells and edges are definitely supported.

FEniCS supports this by using the dx measure for the integral over the triangles, and the dS measure for the integral over the interior edges (facets). You can also use ds (note the small s) for integrals over the domain boundary edges.

There is a demo on the website which uses this functionality at:

http://fenicsproject.org/documentation/dolfin/dev/python/demo/pde/biharmonic/python/documentation.html

Revision history for this message
Claudia Totzeck (totzeck) said : 		#2

Does that mean, that FEniCS treats integrals over the whole domain always as sum of integrals over the triangles? That was not clear to me until now. So thanks for your answer.

In the meantime I found out how to define FunctionSpaces on facets. But I have a new problem now, I want do compute the product of a facet-Function l, a Raviart-Thomas Function t and a facet normal n in this way:

l*dot(t,n)*dS
or
l*inner(t,n)*dS

If I use the dot notation I get the error: " Dot product requires non-scalar arguments, got arguments with ranks 0 and 1."
If I use inner I get : "Shape dismatch."

If I use for debugging grad(t) instead of t, the code works... So I think FEniCS treats a Raviart-Thomas Function not as a vector but a scalar. How can I change this.... I need to use the Function t itset not its gradient...

Thank you!

Revision history for this message
Marie Rognes (meg-simula) said : 		#3

Please provide a minimal, _runnable_ example when asking a specific question: we cannot reproduce your problem without it. So, please provide the definitions of l, t, n etc.

Revision history for this message
Claudia Totzeck (totzeck) said : 		#4

Okay, sorry... Here is the Code:

V = Expression('x[0]*x[0]+x[1]*x[1]')

mesh = UnitSquare(20,20)
DG = FunctionSpace(mesh, "DG", 0)
RT = FunctionSpace(mesh, "RT", 1)
F = FunctionSpace(mesh,'DG',0,'facet')
W = MixedFunctionSpace([RT,DG,F])

#variational form
s = Function(W)
t = TestFunction(W)
d = TrialFunction(W)
n = FacetNormal(mesh)

a = inner(s[0],t[0])*dx - s[2]*inner(t[0],n)*dS
# a = inner(s[0],t[0])*dx - s[2]*dot(t[0],n)*dS

Revision history for this message
Best Marie Rognes (meg-simula) said : 		#5

Try this instead:

--
from dolfin import *

mesh = UnitSquare(20,20)
DG = FunctionSpace(mesh, "DG", 0)
RT = FunctionSpace(mesh, "RT", 1)
F = FunctionSpace(mesh,'DG',0,'facet')
W = MixedFunctionSpace([RT,DG,F])

#variational form
s = Function(W)

# Split s into its 3 components (s0 will be vector-valued)
(s0, s1, s2) = split(s)

# Define separate test and trial functions corresponding to each
# component (often easier to work with)
(t0, t1, t2) = TestFunctions(W)
(d0, d1, d2) = TrialFunctions(W)

n = FacetNormal(mesh)

a = inner(s0, t0)*dx - s2*inner(t0, n)*dS
--

A mixed space will be completely unnested when you refer to the components directly, so when you referred to s[0], you referred to the first component of a RT field. When you use split, you get the "top hierarchy".

Revision history for this message
Claudia Totzeck (totzeck) said : 		#6

Thanks Marie Rognes, that solved my question.