mesh_modification.c

一个用来实现偏微分方程中网格的计算库

				  // New boundary side is Tri 0, side 0
				  new_bndry_elements.push_back(tri0);
				  new_bndry_sides.push_back(0);
				  break;
				}
			      case 1:
				{
				  // New boundary side is Tri 0, side 1
				  new_bndry_elements.push_back(tri0);
				  new_bndry_sides.push_back(1);
				  break;
				}
			      case 2:
				{
				  // New boundary side is Tri 1, side 1
				  new_bndry_elements.push_back(tri1);
				  new_bndry_sides.push_back(1);
				  break;
				}
			      case 3:
				{
				  // New boundary side is Tri 1, side 2
				  new_bndry_elements.push_back(tri1);
				  new_bndry_sides.push_back(2);
				  break;
				}

			      default:
				{
				  std::cerr << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
				  libmesh_error();
				}
			      }
			  }

			else // edge_swap==true
			  {
			    switch (sn)
			      {
			      case 0:
				{
				  // New boundary side is Tri 0, side 0
				  new_bndry_elements.push_back(tri0);
				  new_bndry_sides.push_back(0);
				  break;
				}
			      case 1:
				{
				  // New boundary side is Tri 1, side 0
				  new_bndry_elements.push_back(tri1);
				  new_bndry_sides.push_back(0);
				  break;
				}
			      case 2:
				{
				  // New boundary side is Tri 1, side 1
				  new_bndry_elements.push_back(tri1);
				  new_bndry_sides.push_back(1);
				  break;
				}
			      case 3:
				{
				  // New boundary side is Tri 0, side 2
				  new_bndry_elements.push_back(tri0);
				  new_bndry_sides.push_back(2);
				  break;
				}

			      default:
				{
				  std::cerr << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
				  libmesh_error();
				}
			      }
			  } // end edge_swap==true
		      } // end if (b_id != BoundaryInfo::invalid_id)
		  } // end for loop over sides

		// Remove the original element from the BoundaryInfo structure.
		mesh.boundary_info->remove(*el);

	      } // end if (mesh_has_boundary_data)
	    
	    // Add the newly-created triangles to the temporary vector of new elements.
	    new_elements.push_back(tri0);
	    new_elements.push_back(tri1);

	    // Delete the original element
	    mesh.delete_elem(*el);
	  } // end if (split_elem)
      } // End for loop over elements
  } 

  
    // Now, iterate over the new elements vector, and add them each to
    // the Mesh.
  {
    std::vector<Elem*>::iterator el        = new_elements.begin();
    const std::vector<Elem*>::iterator end = new_elements.end();
    for (; el != end; ++el)
      {
	mesh.add_elem(*el);
      }
  }

  
  if (mesh_has_boundary_data)
    {
      // By this time, we should have removed all of the original boundary sides
      // - except on a hybrid mesh, where we can't "start from a blank slate"! - RHS
      // libmesh_assert (mesh.boundary_info->n_boundary_conds()==0);

      // Clear the boundary info, to be sure and start from a blank slate.
      // mesh.boundary_info->clear();

      // If the old mesh had boundary data, the new mesh better have some.
      libmesh_assert (new_bndry_elements.size() > 0);

      // We should also be sure that the lengths of the new boundary data vectors
      // are all the same.
      libmesh_assert (new_bndry_elements.size() == new_bndry_sides.size());
      libmesh_assert (new_bndry_sides.size()    == new_bndry_ids.size());

      // Add the new boundary info to the mesh
      for (unsigned int s=0; s<new_bndry_elements.size(); ++s)
	mesh.boundary_info->add_side(new_bndry_elements[s],
				     new_bndry_sides[s],
				     new_bndry_ids[s]);
    }
  

  // Prepare the newly created mesh for use.
  mesh.prepare_for_use();

  // Let the new_elements and new_bndry_elements vectors go out of scope.
}


void MeshTools::Modification::smooth (MeshBase& mesh,
                                      const unsigned int n_iterations,
                                      const Real power)
{
  /**
   * This implementation assumes every element "side" has only 2 nodes.
   */
  libmesh_assert (mesh.mesh_dimension() == 2);
  
  /*
   * find the boundary nodes
   */
  std::vector<bool>  on_boundary;
  MeshTools::find_boundary_nodes(mesh, on_boundary);

  for (unsigned int iter=0; iter<n_iterations; iter++)

    {
      /*
       * loop over the mesh refinement level
       */
      unsigned int n_levels = MeshTools::n_levels(mesh);
      for (unsigned int refinement_level=0; refinement_level != n_levels;
           refinement_level++)
        {
          // initialize the storage (have to do it on every level to get empty vectors
          std::vector<Point> new_positions;
          std::vector<Real>   weight;
          new_positions.resize(mesh.n_nodes());
          weight.resize(mesh.n_nodes());

          {
            /*
             * Loop over the elements to calculate new node positions
             */
            MeshBase::element_iterator       el  = mesh.level_elements_begin(refinement_level);
            const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level); 
	
            for (; el != end; ++el)
              {
                /*
                 * Constant handle for the element
                 */
                const Elem* elem = *el;
          
                /*
                 * We relax all nodes on level 0 first 
                 * If the element is refined (level > 0), we interpolate the
                 * parents nodes with help of the embedding matrix
                 */
                if (refinement_level == 0)
                  {
                    for (unsigned int s=0; s<elem->n_neighbors(); s++)
                      {
                        /*
                         * Only operate on sides which are on the
                         * boundary or for which the current element's
                         * id is greater than its neighbor's.
                         * Sides get only built once.
                         */
                        if ((elem->neighbor(s) != NULL) &&
                            (elem->id() > elem->neighbor(s)->id()) ) 
                          {
                            AutoPtr<Elem> side(elem->build_side(s));

                            Node* node0 = side->get_node(0);
                            Node* node1 = side->get_node(1);

                            Real node_weight = 1.;
                            // calculate the weight of the nodes
                            if (power > 0)
                              {
                                Point diff = (*node0)-(*node1);
                                node_weight = std::pow( diff.size(), power );
                              }

                            const unsigned int id0 = node0->id(), id1 = node1->id();
                            new_positions[id0].add_scaled( *node1, node_weight );
                            new_positions[id1].add_scaled( *node0, node_weight );
                            weight[id0] += node_weight;
                            weight[id1] += node_weight;
                          }
                      } // element neighbor loop
                  } 
#ifdef ENABLE_AMR
                else   // refinement_level > 0
                  {
                    /*
                     * Find the positions of the hanging nodes of refined elements.
                     * We do this by calculating their position based on the parent
                     * (one level less refined) element, and the embedding matrix
                     */

                    const Elem* parent = elem->parent();

                    /*
                     * find out which child I am
                     */
                    for (unsigned int c=0; c < parent->n_children(); c++)
                      {
                        if (parent->child(c) == elem)
                          {
                            /*
                             *loop over the childs (that is, the current elements) nodes
                             */
                            for (unsigned int nc=0; nc < elem->n_nodes(); nc++)
                              {
                                /*
                                 * the new position of the node
                                 */
                                Point point;
                                for (unsigned int n=0; n<parent->n_nodes(); n++)
                                  {
                                    /*
                                     * The value from the embedding matrix
                                     */
                                    const float em_val = parent->embedding_matrix(c,nc,n);
	      
                                    if (em_val != 0.)
                                      point.add_scaled (parent->point(n), em_val);
                                  }

                                const unsigned int id = elem->get_node(nc)->id();
                                new_positions[id] = point;
                                weight[id] = 1.;
                              }
                    
                          } // if parent->child == elem
                      } // for parent->n_children
                  } // if element refinement_level
#endif // #ifdef ENABLE_AMR

              } // element loop
          
            /*
             * finally reposition the vertex nodes
             */
            for (unsigned int nid=0; nid<mesh.n_nodes(); ++nid)
              if (!on_boundary[nid] && weight[nid] > 0.)
                mesh.node(nid) = new_positions[nid]/weight[nid];
          }

          {
            /*
             * Now handle the additional second_order nodes by calculating
             * their position based on the vertex postitions
             * we do a second loop over the level elements
             */
            MeshBase::element_iterator       el  = mesh.level_elements_begin(refinement_level);
            const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level); 

            for (; el != end; ++el)
              {
                /*
                 * Constant handle for the element
                 */
                const Elem* elem = *el;
                const unsigned int son_begin = elem->n_vertices();
                const unsigned int son_end   = elem->n_nodes();
                for (unsigned int n=son_begin; n<son_end; n++)
                  {
                    const unsigned int n_adjacent_vertices =
                      elem->n_second_order_adjacent_vertices(n);

                    Point point; 
                    for (unsigned int v=0; v<n_adjacent_vertices; v++)
                      point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));

                    const unsigned int id = elem->get_node(n)->id();
                    mesh.node(id) = point/n_adjacent_vertices;
                  }
              }
          }
      
        } // refinement_level loop

    } // end iteration
}



#ifdef ENABLE_AMR
void MeshTools::Modification::flatten(MeshBase& mesh)
{
  // Algorithm:
  // .) For each active element in the mesh: construct a 
  //    copy which is the same in every way *except* it is
  //    a level 0 element.  Store the pointers to these in
  //    a separate vector. Save any boundary information as well.
  //    Delete the active element from the mesh.
  // .) Loop over all (remaining) elements in the mesh, delete them.
  // .) Add the level-0 copies back to the mesh
  
  // Temporary storage for new element pointers
  std::vector<Elem*> new_elements;

  // BoundaryInfo Storage for element ids, sides, and BC ids
  std::vector<Elem*>              saved_boundary_elements;
  std::vector<short int>          saved_bc_ids;
  std::vector<unsigned short int> saved_bc_sides;

  // Reserve a reasonable amt. of space for each
  new_elements.reserve(mesh.n_active_elem());
  saved_boundary_elements.reserve(mesh.boundary_info->n_boundary_conds());
  saved_bc_ids.reserve(mesh.boundary_info->n_boundary_conds());
  saved_bc_sides.reserve(mesh.boundary_info->n_boundary_conds());
  {
    MeshBase::element_iterator       it  = mesh.active_elements_begin();
    const MeshBase::element_iterator end = mesh.active_elements_end(); 

    for (; it != end; ++it)
      {
	Elem* elem = *it;

	// Make a new element of the same type
	Elem* copy = Elem::build(elem->type()).release();

	// Set node pointers (they point to nodes in the original mesh)
	// The copy's element ID will be set upon adding it to the mesh
	for(unsigned int n=0; n<elem->n_nodes(); n++)
	  copy->set_node(n) = elem->get_node(n);

	// This element could have boundary info as well.  We need
	// to save the (elem, side, bc_id) triples
	for (unsigned int s=0; s<elem->n_sides(); s++)
	    if (elem->neighbor(s) == NULL)
	      {
		short int bc_id = mesh.boundary_info->boundary_id (elem,s);

		if (bc_id != BoundaryInfo::invalid_id)
		  {
		    saved_boundary_elements.push_back(copy);
		    saved_bc_ids.push_back(bc_id);
		    saved_bc_sides.push_back(s);
		  }
	      }

	
	// We're done with this element
	mesh.delete_elem(elem);

	// But save the copy
	new_elements.push_back(copy);
      }
    
    // Make sure we saved the same number of boundary conditions
    // in each vector.
    libmesh_assert (saved_boundary_elements.size() == saved_bc_ids.size());
    libmesh_assert (saved_bc_ids.size()            == saved_bc_sides.size());
  }

  
  // Loop again, delete any remaining elements
  {
    MeshBase::element_iterator       it  = mesh.elements_begin();
    const MeshBase::element_iterator end = mesh.elements_end(); 

    for (; it != end; ++it)
      mesh.delete_elem( *it );
  }

  
  // Add the copied (now level-0) elements back to the mesh
  {
    for (std::vector<Elem*>::iterator it = new_elements.begin();
	 it != new_elements.end();
	 ++it)
      mesh.add_elem(*it);
  }

  // Finally, also add back the saved boundary information
  for (unsigned int e=0; e<saved_boundary_elements.size(); ++e)
    mesh.boundary_info->add_side(saved_boundary_elements[e],
				 saved_bc_sides[e],
				 saved_bc_ids[e]);

  // Trim unused and renumber nodes and elements
  mesh.prepare_for_use();
}
#endif // #ifdef ENABLE_AMR