Subject: [UFC-dev] position of nodes

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 23 Mar 2007 17:24:35 +0100

Is there a function that returns the positions of nodes with respect to 
parent domain? For higher-order elements, I'd like to compute the global 
coordinates of nodes on a facet so that I can determine whether or not 
Dirichlet boundary conditions need to be applied to the node.

Garth

---

From: Anders Logg <[email-masked]>
Date: Fri, 23 Mar 2007 20:05:26 +0100

> Dirichlet boundary conditions need to be applied to the node.

This should not be necessary. The function tabulate_facet_dofs() in
ufc::dof_map() will tabulate the global dofs of a cell that are on a
given facet. The function num_facet_dofs() returns how many facet dofs
there are. So we no longer need to compute the physical coordinates of
the dofs (which only works for point evaluation anyway) and check if
they are on the boundary.

FFC currently generates the following code for tabulate_facet_dofs:

  virtual void tabulate_facet_dofs(unsigned int* dofs,
                                   const ufc::mesh& m,
                                   const ufc::cell& c,
                                   unsigned int facet) const
  {
    // Not implemented
  }

This is no big deal to fix, but it's one of the few functions that
remain to implement (generate).

Also note that Dirichlet boundary condition will mean different things
for different elements. For Lagrange elements, it will mean fixing the
value of a function at the boundary, while for BDM elements it will
mean fixing the normal component at the boundary.

/Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 23 Mar 2007 20:16:09 +0100

> they are on the boundary.

I don't see how this will work if the boundary conditions are provided 
as a function of (x,y,z)?

Garth

> _______________________________________________
> UFC-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufc-dev
>

---

From: Anders Logg <[email-masked]>
Date: Fri, 23 Mar 2007 20:24:47 +0100

> as a function of (x,y,z)?

It will work by evaluating the corresponding dofs (linear functionals)
on the function that defines the boundary condition. This can be done
by calling evaluate_dof() in ufc::finite_element.

So for Lagrange elements, evaluate_dof(i, f, c) will return f(x_i)
where x_i is the i:th node on the cell c.

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Sat, 24 Mar 2007 13:53:41 +0100

> where x_i is the i:th node on the cell c.

But we might still have a problem: boundary conditions aren't
necessarily tied to a facet, they can be tied to f.ex. a vertex as
well. We have no way to set Dirichlet conditions in f.ex. a given
single vertex. Also, when two different conditions apply to
neighbouring boundary facets, the user (on the DOLFIN/PyCC side) must
be free to choose which one applies to the dofs shared by these two
facets. One solution to this without changing UFC would be to let the
user define which BC takes precedence in the case of shared dofs, but
I'm not quite sure if this is flexible enough.

Maybe we should add a few other functions, these might be useful for
other things as well:

class dof_map
{
...

    // Return true if get_coordinate is available, that is
    // if the dofs have associated coordinates
    virtual bool has_coordinates() const  = 0;

    // Return the coordinates associated with dof i on cell c
    virtual void get_coordinates(double *coords, unsigned int i, const
cell &c) = 0;

};

--
---
Martin Sandve Alnæs
PhD student, Cardiac Computations
Simula Research Laboratory, Norway

---

From: Anders Logg <[email-masked]>
Date: Mon, 26 Mar 2007 13:46:52 +0200

>     virtual void get_coordinates(double *coords, unsigned int i, const

I'm not sure we need these. Wouldn't we solve the problem if a
boundary condition must always be specified by a Function (the same
thing that goes into the w array). Then we can have the same problem
(and solution) to both coefficients (like right-hand sides where
someone might also want to do crazy things like specifying the function
at a single vertex) and boundary conditions...

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 26 Mar 2007 14:34:14 +0200

> at a single vertex) and boundary conditions...

I don't see how that solves the problem of actually finding the dof
which is located at a single vertex?

martin

---

From: Anders Logg <[email-masked]>
Date: Mon, 26 Mar 2007 14:51:13 +0200

> which is located at a single vertex?

Fixing the boundary condition at a single vertex does not come natural
in a finite element method. What does it mean mathematically?

Anyway, wouldn't you need to have a mapping from vertices to dofs to
do this? I would suggest that if you want to do something very special
like fixing the values at a single vertex, then you are free to do so
but you need to do something special to do it. For example, you could
read an appendix in the FFC manual (which is not there yet) and see
that the vertex dofs are always ordered first in the dof_map generated
by FFC, so if you want to set a Dirichlet condition at vertex i, then
row i in the matrix and vector is where to set the Dirichlet
condition.

/Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 26 Mar 2007 15:17:49 +0200

Anders Logg wrote:
> in a finite element method. 

Perhaps not, but it is done very often in practice.

What does it mean mathematically?

Usually a Delta function.

Garth

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 26 Mar 2007 15:22:52 +0200

2007/3/26, Anders Logg <[email-masked]>:
>I would suggest that if you want to do something very special

I guess we can't cover everything in a general interface, so for rare
operations this is inevitable. But I think getting the coordinates for
each dof is a very useful operation with many practical uses, which is
trivial to implement as well.

martin

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 26 Mar 2007 15:57:51 +0200

> at a single vertex) and boundary conditions...

Now that I think more about it, we have an application where it's very 
important to know the location of different dofs.

Garth

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 26 Mar 2007 16:32:05 +0200

> important to know the location of different dofs.

Don't know if you're thinking about quadrature?? but I do need the coordinates
of different dofs in order to compute the Jacobian for higher order elements. So
it would be nice if this could be done with a simple dof_map.

Kristian

---

From: Anders Logg <[email-masked]>
Date: Tue, 27 Mar 2007 10:51:04 +0200

On Mon, Mar 26, 2007 at 04:32:05PM +0200, Kristian Oelgaard wrote:
> Don't know if you're thinking about quadrature?? but I do need the
> coordinates of different dofs in order to compute the Jacobian for
> higher order elements. So it would be nice if this could be done
> with a simple dof_map.
> 
> Kristian

Since everyone wants to know where the points are, maybe we should
include it in the interface. :-) Let's add something like

    virtual void tabulate_points(double** coordinates,
                                 const mesh& m,
                                 const cell& c) const = 0;

to the UFC interface. But at the same time, we should try to avoid
using it in DOLFIN if possible.

Maybe we also need something like has_points(), or we could always try
to compute something sensible in tabulate_points() even for elements
that don't have dofs defined by point evaluation. For example, if a
dof is defined by a line integral over an edge, then we could return
the midpoint of the edge in tabulate_points.

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Fri, 30 Mar 2007 15:08:13 +0200

> the midpoint of the edge in tabulate_points.

Ok, and if there is no relation to geometric entities of the cell, the
midpoint of the cell is a good default.

I'll add this function, as well as the sub_dof_map stuff mentioned in
another mail.

martin

---

From: Anders Logg <[email-masked]>
Date: Mon, 2 Apr 2007 13:37:19 +0200

> I'll add this function, as well as the sub_dof_map stuff mentioned in

Did you add tabulate_points? I can't find it.

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 2 Apr 2007 13:52:45 +0200

> Did you add tabulate_points? I can't find it.

It's called tabulate_coordinates. At the end of class dof_map.

martin

---

From: Anders Logg <[email-masked]>
Date: Mon, 2 Apr 2007 13:54:27 +0200

> It's called tabulate_coordinates. At the end of class dof_map.

ok, very good.

/Anders