On Mon, 6 Dec 2004, John Peterson wrote:
> Mark Blome writes:
> >
> > I am currently developing a finite element geoelectric forward
> > solver and I am very pleased with the possibilities libmesh gives
> > me for doing this. For the equation system assembly function I
> > ran across a problem for which I cannot really find a solution.
> > For the geolelectric problem I have to introduce a source term
> > associated to a node on the surface of the input mesh (for the
> > current electrode). The node is correctly marked with a boundary
> > marker and I can read it in without a problem (its a tetgen mesh).
> > However I cannot find a way to set the source term in the right
> > hand side correctly because I can not find out which one of the
> > quadrature points is the one that corresponds to the node that I
> > marked in the mesh. Am I right that it is not even guaranteed that
> > one of the quadrature points will lie on the node where my
> > electrode is attached to? Would that mean that I have to
> > interpolate the source term from the node position to the nearest
> > quadrature point to get the correct right hand side value? Is
> > there a way to force one of the quadrature points to be on the
> > node? Or are my thoughts just going in the wrong direction ?
> > Thanks for any help on this,
> You are correct that the Gaussian quadrature points do not lie on the
> nodes. There are the trapezoidal and Simpson quadrature rules whose points
> however do lie at the nodes. If the value of the source term is availabe
> at every node, you can use the finite element basis functions to interpolate
> a value at the quadrature point. I assume that it makes sense for your
> problem for the source term at a given point to be a weighted average of
> nearby points...
It sounds as if this source term is a Dirac delta functional, not an
integrable function. If that's the case, then unless the functional
happens to fall on a quadrature point it can't be exactly evaluated in
the quadrature loop at all.
I'm assuming the term to be integrated looks like the integral of
k * delta(p) * v_i, and assuming you're using continuous elements.
If the point load is on the interior of an element, you can simply use
the FE<Dim,T>::shape() function to get the value of test function v_i
at the point load location p, and multiply this by the load magnitude
k.
If the point load is on the edge of an element but on the interior of
the domain, then you'll also need to make sure you only add the load
on one element, not once for each adjoining element.
If the point load is on the boundary of your domain, then the right
thing to do depends on your problem physics and on how the load is
defined. In the worst case you may need to figure out what fraction
of the angle (out of 2 pi radians or 4 pi steradians) at that boundary
point is inside your domain, and only add that fraction of the load
magnitude.

Roy Stogner
