--- Begin Message ---
------------------------------------------------------------
revno: 5765
committer: Marie E. Rognes <meg@xxxxxxxxx>
branch nick: dolfin
timestamp: Tue 2011-03-15 20:07:33 +0100
message:
Add barycenter based alternative for finding associated cells.
modified:
dolfin/function/Function.cpp
dolfin/function/Function.h
dolfin/mesh/Mesh.cpp
dolfin/mesh/Mesh.h
--
lp:~dolfin-core/dolfin/rognes
https://code.launchpad.net/~dolfin-core/dolfin/rognes
Your team DOLFIN Core Team is subscribed to branch lp:~dolfin-core/dolfin/rognes.
To unsubscribe from this branch go to https://code.launchpad.net/~dolfin-core/dolfin/rognes/+edit-subscription
=== modified file 'dolfin/function/Function.cpp'
--- dolfin/function/Function.cpp 2011-02-24 17:44:26 +0000
+++ dolfin/function/Function.cpp 2011-03-15 19:07:33 +0000
@@ -6,7 +6,7 @@
// Modified by Andre Massing, 2009.
//
// First added: 2003-11-28
-// Last changed: 2011-02-21
+// Last changed: 2011-03-15
#include <algorithm>
#include <map>
@@ -377,14 +377,52 @@
{
assert(_function_space);
- // Check if UFC cell comes from mesh, otherwise redirect to point-based evaluation
+ // Check if UFC cell comes from mesh, otherwise redirect to
+ // evaluate on non-matching cell
if (ufc_cell.mesh_identifier == (int) _function_space->mesh().id())
{
const Cell cell(_function_space->mesh(), ufc_cell.index);
eval(values, x, cell, ufc_cell);
}
else
- eval(values, x);
+ non_matching_eval(values, x, ufc_cell);
+}
+//-----------------------------------------------------------------------------
+void Function::non_matching_eval(Array<double>& values,
+ const Array<double>& x,
+ const ufc::cell& ufc_cell) const
+{
+ assert(_function_space);
+
+ const double* _x = x.data().get();
+ const Point point(_function_space->mesh().geometry().dim(), _x);
+
+ // Alternative 1: Find cell that point (x) intersects
+ int id = _function_space->mesh().any_intersected_entity(point);
+
+ if (id == -1 && !allow_extrapolation)
+ error("Unable to evaluate function at given point (not inside domain). Set parameter 'allow_extrapolation' to true to allow extrapolation'.");
+
+ // Alternative 2: Compute cell that contains barycenter of ufc_cell
+ // NB: This is slightly heuristic, but should work well for
+ // evaluation of points on refined meshes
+ if (id == -1 && allow_extrapolation)
+ id = _function_space->mesh().supporting_cell(ufc_cell);
+
+ // Alternative 3: Compute closest cell to point (x)
+ if (id == -1 && allow_extrapolation)
+ id = _function_space->mesh().closest_cell(point);
+
+ // Throw error if all alternatives failed.
+ if (id == -1)
+ error("Cannot evaluate function at given point. No matching cell found.");
+
+ // Create cell that contains point
+ const Cell cell(_function_space->mesh(), id);
+ const UFCCell new_ufc_cell(cell);
+
+ // Call evaluate function
+ eval(values, x, cell, new_ufc_cell);
}
//-----------------------------------------------------------------------------
void Function::restrict(double* w,
@@ -396,7 +434,8 @@
assert(_function_space);
// Check if we are restricting to an element of this function space
- if (_function_space->has_element(element) && _function_space->has_cell(dolfin_cell))
+ if (_function_space->has_element(element)
+ && _function_space->has_cell(dolfin_cell))
{
// Get dofmap for cell
const GenericDofMap& dofmap = _function_space->dofmap();
=== modified file 'dolfin/function/Function.h'
--- dolfin/function/Function.h 2011-01-31 17:42:02 +0000
+++ dolfin/function/Function.h 2011-03-15 19:07:33 +0000
@@ -7,7 +7,7 @@
// Modified by Andre Massing, 2009.
//
// First added: 2003-11-28
-// Last changed: 2011-01-31
+// Last changed: 2011-03-15
#ifndef __FUNCTION_H
#define __FUNCTION_H
@@ -140,7 +140,12 @@
virtual uint value_dimension(uint i) const;
/// Evaluate function for given data
- virtual void eval(Array<double>& values, const Array<double>& x, const ufc::cell& cell) const;
+ virtual void eval(Array<double>& values, const Array<double>& x,
+ const ufc::cell& cell) const;
+
+ /// Evaluate function for given data
+ void non_matching_eval(Array<double>& values, const Array<double>& x,
+ const ufc::cell& ufc_cell) const;
/// Restrict function to local cell (compute expansion coefficients w)
virtual void restrict(double* w,
=== modified file 'dolfin/mesh/Mesh.cpp'
--- dolfin/mesh/Mesh.cpp 2011-03-12 21:50:12 +0000
+++ dolfin/mesh/Mesh.cpp 2011-03-15 19:07:33 +0000
@@ -8,7 +8,7 @@
// Modified by Andre Massing, 2009-2010.
//
// First added: 2006-05-09
-// Last changed: 2011-02-07
+// Last changed: 2011-03-15
#include <dolfin/ale/ALE.h>
#include <dolfin/common/Timer.h>
@@ -370,6 +370,24 @@
return _intersection_operator.closest_point_and_cell(point);
}
//-----------------------------------------------------------------------------
+dolfin::uint Mesh::supporting_cell(const ufc::cell& cell) const
+{
+ // Extract vertices of cell from some (other) mesh
+ const double * const * vertices = cell.coordinates;
+
+ // Compute barycenter
+ Point barycenter;
+ for (uint i = 0; i <= _topology.dim(); i++)
+ {
+ Point vertex(_topology.dim(), vertices[i]);
+ barycenter += vertex;
+ }
+ barycenter /= (_topology.dim() + 1);
+
+ // Find cell intersecting this barycenter
+ return any_intersected_entity(barycenter);
+}
+//-----------------------------------------------------------------------------
IntersectionOperator& Mesh::intersection_operator()
{
return _intersection_operator;
=== modified file 'dolfin/mesh/Mesh.h'
--- dolfin/mesh/Mesh.h 2011-03-12 21:50:12 +0000
+++ dolfin/mesh/Mesh.h 2011-03-15 19:07:33 +0000
@@ -9,7 +9,7 @@
// Modified by Andre Massing, 2009-2010.
//
// First added: 2006-05-08
-// Last changed: 2011-03-10
+// Last changed: 2011-03-15
#ifndef __MESH_H
#define __MESH_H
@@ -17,6 +17,7 @@
#include <string>
#include <utility>
#include <boost/scoped_ptr.hpp>
+#include <ufc.h>
#include <dolfin/ale/ALEType.h>
#include <dolfin/common/types.h>
@@ -549,6 +550,20 @@
/// 1
dolfin::uint closest_cell(const Point& point) const;
+ /// Computes the index of the cell in the mesh which contained the
+ /// barycenter of the given ufc_cell
+ ///
+ /// *Arguments*
+ /// cell (_ufc::cell_)
+ /// A ufc::cell object.
+ ///
+ /// *Returns*
+ /// uint
+ /// The index of the cell in the mesh which supports the
+ /// barycenter
+ ///
+ dolfin::uint supporting_cell(const ufc::cell& cell) const;
+
/// Computes the point inside the mesh and the corresponding cell
/// index which are closest to the point query.
///
--- End Message ---
