--- Begin Message ---
------------------------------------------------------------
revno: 4535
committer: Anders Logg <logg@xxxxxxxxx>
branch nick: dolfin-dev
timestamp: Wed 2010-02-17 15:34:42 +0100
message:
Add support for specifying facet orientation in assembly over interior facets
modified:
ChangeLog
dolfin/fem/Assembler.cpp
dolfin/mesh/Facet.h
dolfin/mesh/MeshData.h
--
lp:dolfin
https://code.launchpad.net/~dolfin-core/dolfin/main
Your team DOLFIN Core Team is subscribed to branch lp:dolfin.
To unsubscribe from this branch go to https://code.launchpad.net/~dolfin-core/dolfin/main/+edit-subscription.
=== modified file 'ChangeLog'
--- ChangeLog 2010-02-15 22:12:10 +0000
+++ ChangeLog 2010-02-17 14:34:42 +0000
@@ -1,3 +1,4 @@
+ - Add support for specifying facet orientation in assembly over interior facets
- Allow user to choose which LU package PETScLUSolver uses
- Add computation of intersection between arbitrary mesh entities
- Random access to MeshEntitiyIterators
=== modified file 'dolfin/fem/Assembler.cpp'
--- dolfin/fem/Assembler.cpp 2010-02-16 10:59:43 +0000
+++ dolfin/fem/Assembler.cpp 2010-02-17 14:34:42 +0000
@@ -6,7 +6,7 @@
// Modified by Kent-Andre Mardal, 2008
//
// First added: 2007-01-17
-// Last changed: 2010-02-08
+// Last changed: 2010-02-17
#include <dolfin/log/dolfin_log.h>
#include <dolfin/common/Timer.h>
@@ -272,6 +272,12 @@
mesh.init(mesh.topology().dim() - 1, mesh.topology().dim());
assert(mesh.ordered());
+ // Get interior facet directions (if any)
+ MeshFunction<uint>* facet_orientation = mesh.data().mesh_function("facet orientation");
+ if (facet_orientation && facet_orientation->dim() != mesh.topology().dim() - 1)
+ error("Expecting facet orientation to be defined on facets (not dimension %d).",
+ facet_orientation);
+
// Assemble over interior facets (the facets of the mesh)
Progress p(AssemblerTools::progress_message(A.rank(), "interior facets"), mesh.num_facets());
for (FacetIterator facet(mesh); !facet.end(); ++facet)
@@ -297,7 +303,7 @@
if (!integral) continue;
// Get cells incident with facet
- std::pair<const Cell, const Cell> cells = facet->adjacent_cells();
+ std::pair<const Cell, const Cell> cells = facet->adjacent_cells(facet_orientation);
const Cell& cell0 = cells.first;
const Cell& cell1 = cells.second;
=== modified file 'dolfin/mesh/Facet.h'
--- dolfin/mesh/Facet.h 2009-11-04 10:52:55 +0000
+++ dolfin/mesh/Facet.h 2010-02-17 14:34:42 +0000
@@ -4,7 +4,7 @@
// Modified by Garth N. Wells, 2009.
//
// First added: 2006-06-02
-// Last changed: 2009-11-04
+// Last changed: 2010-02-17
#ifndef __FACET_H
#define __FACET_H
@@ -12,6 +12,7 @@
#include "Cell.h"
#include "MeshEntity.h"
#include "MeshEntityIterator.h"
+#include "MeshFunction.h"
namespace dolfin
{
@@ -29,7 +30,7 @@
~Facet() {}
// FIXME: This function should take care of facet 'ownership' when a mesh
- // is distributed across processes
+ // is distributed across processes
/// Determine whether or not facet is an interior facet. This is 'relative'
/// to the given partition of the mesh if the mesh is distributed
bool interior() const
@@ -43,13 +44,33 @@
}
// FIXME: This function should take care of the case where adjacent cells
- // live on different processes
- /// Return adjacent cells
- std::pair<const Cell, const Cell> adjacent_cells() const
- {
+ // FIXME: live on different processes
+
+ /// Return adjacent cells. An optional argument that lists for
+ /// each facet the index of the first cell may be given to specify
+ /// the ordering of the two cells. If not specified, the ordering
+ /// will depend on the (arbitrary) ordering of the mesh
+ /// connectivity.
+ std::pair<const Cell, const Cell> adjacent_cells(MeshFunction<uint>* facet_orientation=0) const
+ {
assert(num_entities(dim() + 1) == 2);
- return std::make_pair(Cell(mesh(), entities(dim() + 1)[0]),
- Cell(mesh(), entities(dim() + 1)[1]));
+
+ // Get cell indices
+ const uint D = dim() + 1;
+ const uint c0 = entities(D)[0];
+ const uint c1 = entities(D)[1];
+
+ // Normal ordering
+ if (!facet_orientation || (*facet_orientation)[*this] == c0)
+ return std::make_pair(Cell(mesh(), c0), Cell(mesh(), c1));
+
+ // Sanity check
+ if ((*facet_orientation)[*this] != c1)
+ error("Illegal facet orientation specified, cell %d is not a neighbor of facet %d.",
+ (*facet_orientation)[*this], index());
+
+ // Opposite ordering
+ return std::make_pair(Cell(mesh(), c1), Cell(mesh(), c0));
}
};
=== modified file 'dolfin/mesh/MeshData.h'
--- dolfin/mesh/MeshData.h 2010-02-17 11:45:22 +0000
+++ dolfin/mesh/MeshData.h 2010-02-17 14:34:42 +0000
@@ -31,6 +31,10 @@
///
/// "exterior facet domains" - MeshFunction<uint> of dimension D - 1
///
+ /// Facet orientation (used for assembly over interior facets)
+ ///
+ /// "facet orientation" - MeshFunction<uint> of dimension D - 1
+ ///
/// Boundary extraction
///
/// "vertex map" - MeshFunction<uint> of dimension 0
--- End Message ---
