dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

Looks good at first glance. I'll take a closer look before I add it.

-- 
Anders


On Thu, Jan 08, 2009 at 02:23:53PM -0700, Bartosz Sawicki wrote:
> 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();

> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev

Attachment: signature.asc
Description: Digital signature