dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #11590
Re: new refinement method
Please find attached patch with new iterative refinement algorithm.
The iteration process is robust, usually less that 10 iterations is
needed to converge. I've tested it with triangular and tetrahedral
meshes up to 100k cells.
This is my first bigger commit into the project. So please read it
carefully. I tried to write fast and efficient code, but maybe someone
has any ideas how to improve it. Waiting for any opinions...
regrds.
BArtek
On 07/01/09 04:15 PM, Bartosz Sawicki wrote:
On 31/12/08 12:51 AM, Anders Logg wrote:
On Tue, Dec 30, 2008 at 04:54:05PM -0700, Bartosz Sawicki wrote:
I'm starting to work on implementation of more advanced algorithm for
mesh refinement. At this moment the simplest version of the
Longest-Edge-Bisection is used inside the Mesh::refine().
I plan to develop some algorithm based on Edge-Bisection with two set of
markers and iterative swapping. There are some papers about this method.
It is thought to be robust and produce good quality mesh.
This sounds good. Our current refinement algorithm sometimes produces
bad quality mesh (very thin triangles) if called many times without
intermediate calls to mesh.smooth().
Today I was stopped by the first problem related with edge global index.
We know that every bisection of triangle produces 3 new edges, and
removes 1 old edge. The question is, if there is method to find what
will be the global index of the new edges. If I understand correctly
edge index can created (mesh.init(1)) only after mesh object is created
and fulfilled by Vertexes and Cells.
Maybe I'm wrong and there is simple way to get edge index even without
calling mesh.init(1)?
Unfortunately no. It would be very good if there were a way to know
global edge numbers. It would also help with the parallel implementation
currently in progress.
But there might be some clever numbering strategy that could work. If
anyone knows of such a numbering strategy (uniquely assigning global
edge numbers based on global vertex numbers) we could implement that
as a post-processing step in TopologyComputation.
One first suggestion would be to numberthe first edge of the first
vertex 0, then the second edge connected to the first vertex etc.
I have finished first working version of new mesh refinement algorithm.
The code is quite dirty now, but I will try to make it more readable for
everyone. I'll send you a patch.
Edge numbering problem was solved by using std:map<Array<int>, int,
cmp>. I've used two vertex number as a no-changing edge identifier.
This works quite fast, but I'm open for new ideas. Maybe someone has
better idea how to build MeshFunction connected with edges, and
transform it into different mesh?
Similar problem would be crucial for implementation of MeshFunctions
transformations connected with mesh refinement. As boundary indicators,
and material properties should be preserved when we change the mesh.
I plan to implement this in the next step.
regrd.
BArtek
------------------------------------------------------------------------
_______________________________________________
DOLFIN-dev mailing list
DOLFIN-dev@xxxxxxxxxx
http://www.fenics.org/mailman/listinfo/dolfin-dev
_______________________________________________
DOLFIN-dev mailing list
DOLFIN-dev@xxxxxxxxxx
http://www.fenics.org/mailman/listinfo/dolfin-dev
diff -r 946e884ef099 dolfin/mesh/LocalMeshRefinement.cpp
--- a/dolfin/mesh/LocalMeshRefinement.cpp Wed Jan 07 23:10:49 2009 +0000
+++ b/dolfin/mesh/LocalMeshRefinement.cpp Thu Jan 08 12:56:06 2009 -0700
@@ -19,6 +19,18 @@
#include "LocalMeshRefinement.h"
using namespace dolfin;
+
+struct cmp2
+{
+ bool operator()(Array<uint> const a, Array<uint> const b)
+ {
+ if(a[0] == b[0]) {
+ return a[1] < b[1];
+ } else {
+ return a[0] < b[0];
+ }
+ }
+};
//-----------------------------------------------------------------------------
void LocalMeshRefinement::refineMeshByEdgeBisection(Mesh& mesh,
@@ -236,9 +248,30 @@
mesh = refined_mesh;
}
//-----------------------------------------------------------------------------
+void LocalMeshRefinement::refineIterativelyByEdgeBisection(Mesh& mesh,
+ MeshFunction<bool>& cell_marker)
+{
+
+ begin("Iterative mesh refinement algorithm");
+
+ MeshFunction<uint> edges(mesh, 1);
+ edges = 0;
+
+ int i = 0;
+ bool more_iterations = true;
+ while( more_iterations )
+ {
+ cout << "=== Iteration : " << i++ << " ===" << endl;
+ more_iterations = iterationOfRefinement(mesh, cell_marker, edges);
+ }
+
+ end();
+
+}
+//-----------------------------------------------------------------------------
void LocalMeshRefinement::bisectEdgeOfSimplexCell(const Cell& cell,
Edge& edge,
- uint& new_vertex,
+ uint new_vertex,
MeshEditor& editor,
uint& current_cell)
{
@@ -272,4 +305,182 @@
}
//-----------------------------------------------------------------------------
+bool LocalMeshRefinement::iterationOfRefinement(Mesh& mesh,
+ MeshFunction<bool>& cell_marker,
+ MeshFunction<uint>& bisected_edges
+ )
+{
+ // Map used for edge function transformation
+ std::map<Array<uint>, uint, cmp2> edge_map;
+ std::map<Array<uint>, uint, cmp2>::iterator edge_map_it;
+ const uint num_vertices = mesh.size(0);
+ const uint num_cells = mesh.size(mesh.topology().dim());
+ const CellType& cell_type = mesh.type();
+
+ Mesh refined_mesh;
+ MeshEditor editor;
+ editor.open(refined_mesh, cell_type.cellType(),
+ mesh.topology().dim(), mesh.geometry().dim());
+
+ // Generate cell - edge connectivity if not generated
+ mesh.init(mesh.topology().dim(), 1);
+ // Generate edge - vertex connectivity if not generated
+ mesh.init(1, 0);
+
+ uint num_new_vertices = 0;
+ uint num_new_cells = 0;
+ uint longest_edge_index = 0;
+ double lmax, l;
+
+ // Temporary edge marker function
+ MeshFunction<bool> found_edge(mesh, 1);
+ found_edge = false;
+
+ // Calculate number of new vertices
+ for (CellIterator c(mesh); !c.end(); ++c)
+ {
+ if( cell_marker.get(c->index()) )
+ {
+ // Find longest edge of cell c
+ lmax = 0.0;
+ for(EdgeIterator e(*c); !e.end(); ++e)
+ {
+ l = e->length();
+ if ( lmax < l )
+ {
+ lmax = l;
+ longest_edge_index = e->index();
+ }
+ }
+ Edge longest_edge(mesh,longest_edge_index);
+ uint middle_vertex = bisected_edges.get(longest_edge_index);
+ if( middle_vertex == 0 && !found_edge.get(longest_edge_index))
+ {
+ found_edge.set(longest_edge_index, true);
+ num_new_vertices++;
+ }
+ num_new_cells++;
+ }
+ }
+
+ cout << "Number of cells in old mesh: " << num_cells << "; to add: " << num_new_cells << endl;
+ cout << "Number of vertices in old mesh: " << num_vertices << "; to add: " << num_new_vertices << endl;
+
+ editor.initVertices(num_vertices + num_new_vertices);
+ editor.initCells(num_cells + num_new_cells);
+
+ //Rewrite old vertices
+ uint current_vertex = 0;
+ for (VertexIterator v(mesh); !v.end(); ++v)
+ editor.addVertex(current_vertex++, v->point());
+
+ //Rewrite old cells
+ uint current_cell = 0;
+ Array<uint> cell_vertices(cell_type.numEntities(0));
+ for (CellIterator c(mesh); !c.end(); ++c)
+ {
+ if( ! cell_marker.get(c->index()) )
+ {
+ uint cv = 0;
+ for (VertexIterator v(*c); !v.end(); ++v)
+ cell_vertices[cv++] = v->index();
+ editor.addCell(current_cell++, cell_vertices);
+ }
+ }
+
+ //Main loop
+ for (CellIterator c(mesh); !c.end(); ++c)
+ {
+ if( cell_marker.get(c->index()) )
+ {
+
+ // Find longest edge of cell c
+ lmax = 0.0;
+ for(EdgeIterator e(*c); !e.end(); ++e)
+ {
+ l = e->length();
+ if ( lmax < l )
+ {
+ lmax = l;
+ longest_edge_index = e->index();
+ }
+ }
+
+ Edge longest_edge(mesh,longest_edge_index);
+
+ uint middle_vertex = bisected_edges.get(longest_edge_index);
+ if( middle_vertex == 0)
+ {
+ // Add new vertex
+ editor.addVertex(current_vertex++, longest_edge.midpoint());
+
+ middle_vertex = current_vertex-1;
+
+ // Insert new vertex into the bisected edges mapping
+ Array<uint> ev(2);
+ ev[0] = longest_edge.entities(0)[0];
+ ev[1] = longest_edge.entities(0)[1];
+ edge_map[ev] = middle_vertex;
+
+ bisected_edges.set(longest_edge_index, middle_vertex);
+ }
+
+ // Add two new cells
+ bisectEdgeOfSimplexCell(*c, longest_edge, middle_vertex+1, editor, current_cell);
+
+ }
+ }
+
+ editor.close();
+
+ // Rewrite old bisected edges
+ for (EdgeIterator e(mesh); !e.end(); ++e)
+ {
+ uint eix = bisected_edges.get(e->index());
+ if( eix )
+ {
+ Array<uint> ev(2);
+ ev[0] = e->entities(0)[0];
+ ev[1] = e->entities(0)[1];
+ edge_map[ev] = eix;
+ }
+ }
+
+ // Initialize new cell markers
+ cell_marker.init(refined_mesh, refined_mesh.topology().dim());
+ cell_marker = false;
+
+ // Initialize new bisected edges
+ bisected_edges.init(refined_mesh, 1);
+ bisected_edges = 0;
+
+ Progress p("Calculate cell markers and bisected edges", refined_mesh.numEdges());
+ bool next_iteration = false;
+ for (EdgeIterator e(refined_mesh); !e.end(); ++e)
+ {
+ Array<uint> ev(2);
+ ev[0]= e->entities(0)[0];
+ ev[1]= e->entities(0)[1];
+ if ( (edge_map_it = edge_map.find(ev)) == edge_map.end())
+ {
+ bisected_edges.set(e->index(), 0);
+ }
+ else
+ {
+ bisected_edges.set(e->index(), edge_map_it->second);
+ for(CellIterator c(*e); !c.end(); ++c) {
+ cell_marker.set(c->index(), true);
+ next_iteration = true;
+ }
+ }
+ p++;
+ }
+
+ // Overwrite old mesh
+ mesh = refined_mesh;
+
+ return next_iteration;
+}
+//-----------------------------------------------------------------------------
+
diff -r 946e884ef099 dolfin/mesh/LocalMeshRefinement.h
--- a/dolfin/mesh/LocalMeshRefinement.h Wed Jan 07 23:10:49 2009 +0000
+++ b/dolfin/mesh/LocalMeshRefinement.h Thu Jan 08 12:56:06 2009 -0700
@@ -26,14 +26,21 @@
static void refineMeshByEdgeBisection(Mesh& mesh,
MeshFunction<bool>& cell_marker,
bool refine_boundary = true);
-
+
+ /// Iteratively refine mesh locally by the longest edge bisection
+ static void refineIterativelyByEdgeBisection(Mesh& mesh,
+ MeshFunction<bool>& cell_marker);
private:
/// Bisect edge of simplex cell
static void bisectEdgeOfSimplexCell(const Cell& cell, Edge& edge,
- uint& new_vertex,
+ uint new_vertex,
MeshEditor& editor,
uint& current_cell);
+
+ static bool iterationOfRefinement(Mesh& mesh,
+ MeshFunction<bool>& cell_marker,
+ MeshFunction<uint>& bisected_edges);
};
diff -r 946e884ef099 dolfin/mesh/Mesh.cpp
--- a/dolfin/mesh/Mesh.cpp Wed Jan 07 23:10:49 2009 +0000
+++ b/dolfin/mesh/Mesh.cpp Thu Jan 08 12:56:06 2009 -0700
@@ -201,6 +201,14 @@
_ordered = false;
}
//-----------------------------------------------------------------------------
+void Mesh::refine_iteratively(MeshFunction<bool>& cell_markers)
+{
+ LocalMeshRefinement::refineIterativelyByEdgeBisection(*this, cell_markers);
+
+ // Mesh may not be ordered
+ _ordered = false;
+}
+//-----------------------------------------------------------------------------
void Mesh::coarsen()
{
// FIXME: Move implementation to separate class and just call function here
diff -r 946e884ef099 dolfin/mesh/Mesh.h
--- a/dolfin/mesh/Mesh.h Wed Jan 07 23:10:49 2009 +0000
+++ b/dolfin/mesh/Mesh.h Thu Jan 08 12:56:06 2009 -0700
@@ -146,8 +146,12 @@
/// Refine mesh uniformly
void refine();
- /// Refine mesh according to cells marked for refinement
+ /// Refine mesh according to cells marked for refinement,
void refine(MeshFunction<bool>& cell_markers, bool refine_boundary = true);
+
+ /// Refine mesh according to cells marked for refinement,
+ /// Iterative the longest edge bisection algorithm
+ void refine_iteratively(MeshFunction<bool>& cell_markers);
/// Coarsen mesh uniformly
void coarsen();
Follow ups
References
