I am trying to integrate the product of two Gridfunction over a specific element. I know we can achieve this by setting “element_wise=True” as in " Integrate(gf1*gf2,mesh,element_wise=True)".
However, this easy option computes the integral value over all the elements and returns a vector made of all these values. This doesn’t seem to be very optimal in computational efficiency if I only need the integral value over one specific element.
I am wondering if there is any other convenient way to compute such an integral over only one element using NGSolve.
you can use a BitArray to mark that one element and integrate using the definedonelements flag. The code would look something like
marker = BitArray(mesh.ne) # mesh.ne = number of elements
marker.Clear()
marker.Set(i) # with i the number of the element
Integrate(gf1 * gf2 * dx(defineonelements=marker), mesh)
As far as I know, this is an efficient way to only integrate over the marked elements.
Thanks a lot for your help. I think using “Integrate(gf1 * gf2 * dx(defineonelements=marker), mesh)” doesn’t work for me, while the symbolic integrator/bilinear forms works well as follows.
mesh = Mesh(unit_square.GenerateMesh(maxh=0.5))
marker = BitArray(mesh.ne) # mesh.ne = number of elements
marker.Clear()
marker.Set(4) # with i the number of the element
U = H1(mesh, order=1)
gf = GridFunction(U)
gf.Set(1)
a = BilinearForm(U)
a += U.TrialFunction()*U.TestFunction() * dx(definedonelements=marker)
print(2*a.Energy(gf.vec))
Moreover, I am wondering if there is a similar convenient way that we can integrate over a specific edge by any chance. Is there something like ds(…) specifying the integrator on one edge? Many thanks again.
Strange, for me Integrate(gf1 * gf2 * dx(defineonelements=marker), mesh) works totally fine. You can add element_boundary=True inside the DifferntialSymbol (dx). Then, you integrate over the boundaries of that element (i.e. the facets). I don’t know how to integrate only over one specific facet using a marker (ds integrates over the boundary of the mesh, not the given elements). But you can define a facet indicator function via GridFunction(FESpace("facet", mesh, order=0)) and set the values on all facets to zero besides the ones you want to integrate on. Multiping this indicator with your function and only integrating over the edges of the specific elements should lead to the desired output.
Hi,
it’s not a single Python statement, but you can use an integration rule for that:
from ngsolve import *
mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))
ir = IntegrationRule(TRIG,5)
cf = x
globigl = 0
for el in mesh.Elements(VOL):
trafo = mesh.GetTrafo(el)
igt = sum( ip.weight*trafo(ip).measure*cf(trafo(ip)) for ip in ir )
globigl += igt
print (igt)
print ("global integral", globigl)
Another idea is to compute the element-vector for a finite element element with one shape-function constant 1 (i.e. an L2(order=0) -finite element):
fes = L2(mesh, order=0)
v = fes.TestFunction()
lfi = (cf*v*dx)[0].MakeLFI()
globigl = 0
for el in mesh.Elements(VOL):
trafo = mesh.GetTrafo(el)
fel = fes.GetFE(el)
igt = lfi.CalcElementVector(fel, trafo)[0]
globigl += igt
print (igt)
print ("global integral", globigl)
I am also wondering if there is any convenient way to set up the quadrature rule on edges so that we can compute the integral over one specific edge. Is it possible to have an analogue of your first solution in the case of integrals on edges? Or should I get the coordinates of points on the edge and set the quadrature rule manually?