dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Thu, Feb 11, 2010 at 08:32:17AM -0800, Johan Hake wrote:
> On Thursday 11 February 2010 06:31:18 dbeacham wrote:
> > Question #99875 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/99875
> >
> >     Status: Answered => Solved
> >
> > dbeacham confirmed that the question is solved:
> > Thanks, this had put me on the right track. Just one more thing
> > (hopefully):
> >
> > In the python interface, how do I generate the 'ufc::cell const &'
> > argument used in DofMap.local_dimension? I've tried using
> > V.cell()/V.ufl_element() but they don't work, so I assume they're the
> > incorrect way.
>
> I guess you can't as it is not possible to create an ufc::cell in the Python
> interface (as far as I know...). I am also not sure we want to expose more of
> the UFC interface to PyDOLFIN, and particular not ufc::cell, as it comes with
> some nasty double ** attributes in its public interface.
>
> However, in DofMap many functions come in one DOLFIN and one UFC version, but
> local_dimension is not one of them. Is there a reason for this Anders? It
> should be straight forward to just add it.
>
> Johan

It is actually possible to call some of them. I have added
dolfin::Cell versions for some of UFC functions, for example
tabulate_dofs in DofMap:

  /// Tabulate the local-to-global mapping of dofs on a cell (UFC cell version)
  void tabulate_dofs(uint* dofs, const ufc::cell& ufc_cell, uint cell_index) const;

  /// Tabulate the local-to-global mapping of dofs on a cell (DOLFIN cell version)
  void tabulate_dofs(uint* dofs, const Cell& cell) const;

I don't think it would be a problem to add more of these. It's useful
for testing out algorithms in Python. See for example the attached
script which makes heavy use of low-level UFC functions from Python.
The algorithm in question has now been ported to C++ in the DOLFIN
Extrapolation class (which was something like a factor 1000 faster).

But in your case, are you sure you need to call DofMap::local_dimension? If the
dimension does not vary over the mesh, you can call
DofMap::max_local_dimension.

--
Anders

__author__ = "Anders Logg (logg@xxxxxxxxx)"
__copyright__ = "Copyright (C) 2009 - Anders Logg"
__license__  = "GNU GPL version 3 or any later version"

# Modified by Marie Rognes, 2009

# Last changed: 2010-02-10

__all__ = ["reconstruct", "do_the_boundary_fix"]

from dolfin import FunctionSpace, Function, Cell, cells, facets, info, DOLFIN_EPS
import numpy
import time

from adaptivity.bernstein import Bernstein

def reconstruct(v, W, domain=None):
    "Create reconstruction of v in FunctionSpace W."

    info("Reconstructing")
    start = time.time()

    mesh = v.function_space().mesh()

    # Temporary 'matrix' to be used to hold cell-wise values before
    # smoothing
    dof_values = [[] for i in range(W.dim())]

    # Array to hold local-to-global dof indices
    dofmap = W.dofmap()
    ltg_dofs = numpy.zeros(dofmap.max_local_dimension(), dtype='I')

    # Step 1: Compute local reconstruction and store values
    for cell in cells(mesh):
        w_T = _compute_local_reconstruction(v, cell, W)

        # Place local dofs in appropriate global dof bucket
        dofmap.tabulate_dofs(ltg_dofs, cell)
        for i in range(len(w_T)):
            dof_values[ltg_dofs[i]] += [w_T[i]]


    # Step 2: Construct reconstruction by averaging local reconstructions
    w = Function(W)
    w.vector()[:] = numpy.array([numpy.mean(d) for d in dof_values])

    # Need to initalize facet-to-cell connectivity:
    d = W.mesh().geometry().dim()
    W.mesh().init(d-1, d)

    # Do the boundary fix
    if domain is not None:
        do_the_boundary_fix(w, domain)

    stop = time.time()
    info("Done! Reconstructing took %d s." % (stop - start))

    return w

def do_the_boundary_fix(w, domain):
    print "Fixing boundary, old-school"
    boundary_dofs = _get_boundary_dofs(w.function_space(), domain)
    w.vector()[boundary_dofs] = numpy.zeros((len(boundary_dofs),), dtype='I')

def _get_boundary_dofs(V, domain):

    domain_dofs = []

    # Need to initalize facet-to-cell connectivity:
    d = V.mesh().geometry().dim()
    V.mesh().init(d-1, d)
    num_facets = d+1

    for facet in facets(V.mesh()):
        if facet.num_entities(d) == 1:

            # The facet is on the boundary and has 1 incident cell
            cell = Cell(V.mesh(), facet.entities(d)[0])
            local_facet_index = cell.index(facet)

            # Get local indices for dofs living on facet
            num_facet_dofs =  V.dofmap().num_facet_dofs()
            local_facet_dofs = numpy.zeros((num_facet_dofs,), dtype='I')
            V.dofmap().tabulate_facet_dofs(local_facet_dofs, local_facet_index)

            # Get the global dof indices of all dofs on cell but
            # choose only those on the facet
            n = V.dofmap().max_local_dimension()
            all_dofs = numpy.zeros((n,), dtype='I')
            V.dofmap().tabulate_dofs(all_dofs, cell)
            dofs = all_dofs[local_facet_dofs]

            # Tabulate coordinates of all dofs on cell
            dof_coordinates = numpy.zeros((n, d))
            V.dofmap().tabulate_coordinates(dof_coordinates, cell)

            for (i, x) in enumerate(dof_coordinates[local_facet_dofs]):
                if domain.inside(x, True):
                    domain_dofs += [dofs[i]]

    return numpy.array(list(set(domain_dofs)))


def _compute_local_reconstruction(v, cell, W):

    # FIXME: Only works for point-values, we should really extract the
    # degrees of freedom (as functionals) associated with the patch:
    # L = extract_dofs(cell/patch) ({L_i}_{i=1}^m):

    # Extract dof coordinates and dof values associated with patch
    points, b = _extract_patch_values(v, cell)
    m = len(points)

    # FIXME: Only works for scalars -- we really need basis for the
    # reconstruction space restricted to the cell (but living on the
    # patch): {phi_j}_{j=1}^n:

    # Create Bernstein basis
    bernstein = Bernstein(degree=W.ufl_element().degree())
    dofmap = W.dofmap()
    n = dofmap.max_local_dimension()

    # Construct local linear system A = L_i(phi[j]):
    A = numpy.array([[bernstein(j, points[i]) for j in range(n)]
                     for i in range(m)])

    # Solve linear system for coefficients in Bernstein expansion
    if m > n:
        coefficients, residuals, rank, singular_values = numpy.linalg.lstsq(A, b)
    elif m < n:
        raise RuntimeError, "Linear system for reconstruction under-constrained (%d x %d)." % (m, n)
    else:
        coefficients = numpy.linalg.solve(A, b)

    # FIXME: This part should not be necessary if we use the nodal
    # basis on T directly (instead of Bernstein). Now need to go from
    # Bernstein expansion to nodal basis expansion

    # FIXME: Want local degrees of freedom for W ({M_j}_{j=1}^n)
    d = W.mesh().geometry().dim()
    dof_coordinates = numpy.zeros((n, d))
    dofmap.tabulate_coordinates(dof_coordinates, cell)

    # M = M_j(phi_i)
    M = numpy.array([[bernstein(j, dof_coordinates[i]) for j in range(n)]
                       for i in range(n)])

    # w_T contains the coefficients for the nodal expansion on this cell
    w_T = numpy.dot(M, coefficients)

    return w_T


def _extract_patch_values(v, cell):
    """Extract points and values for v on a patch of elements
    surrounding the given cell."""

    # Note: We assume that the dofs are associated with point
    # evaluation (not moments or normal components etc). It might also
    # only work for scalars at this point. Haven't thought this
    # through in full detail.

    # Get number of dofs on cell and geometric dimension
    dofmap = v.function_space().dofmap()
    mesh = v.function_space().mesh()
    n = dofmap.max_local_dimension()
    d = mesh.geometry().dim()

    # Initialize arrays
    cell_dofs = numpy.zeros(n, dtype='I')
    cell_coordinates = numpy.zeros((n, d))

    # Create list of cell itself and surrounding cells
    cell_indices = [cell.index()] + [c.index() for c in cells(cell)]

    # Get dofs and coordinates
    dofs = []
    coordinates = []
    for cell_index in cell_indices:

        # Create cell
        c = Cell(mesh, cell_index)

        # Get dofs and coordinates
        dofmap.tabulate_dofs(cell_dofs, c)
        dofmap.tabulate_coordinates(cell_coordinates, c)

        # NB: Must copy dofs and coordinates or we will overwrite previous values!

        # Store dofs and coordinates
        for dof in cell_dofs:
            dofs.append(dof.copy())
        for x in cell_coordinates:
            coordinates.append(x.copy())

    # Get values for dofs
    x = v.vector()
    values = [x[dof] for dof in dofs]

    #for i in range(len(values)):
    #    print "  v(%s) = %g" % (", ".join([str(float(c)) for c in coordinates[i]]), values[i])

    return coordinates, values

Attachment: signature.asc
Description: Digital signature