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 ([email protected])" __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
signature.asc
Description: Digital signature
_______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : [email protected] Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
