inf_elem_builder.c

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

    MeshBase::element_iterator       it  = this->_mesh.active_elements_begin();
    const MeshBase::element_iterator end = this->_mesh.active_elements_end();
  
    for(; it != end; ++it)
      {
	Elem* elem = *it;
	
	for (unsigned int s=0; s<elem->n_neighbors(); s++)
	  {
	    // check if elem(e) is on the boundary
	    if (elem->neighbor(s) == NULL)
	      {
		// note that it is safe to use the Elem::side() method, 
		// which gives a non-full-ordered element 
		AutoPtr<Elem> side(elem->build_side(s));

		// bool flags for symmetry detection
		bool sym_side=false;		
		bool on_x_sym=true;
		bool on_y_sym=true;
		bool on_z_sym=true;
		
		
		// Loop over the nodes to check whether they are on the symmetry planes,
		// and therefore sufficient to use a non-full-ordered side element
		for(unsigned int n=0; n<side->n_nodes(); n++)
		  {
		    const Point dist_from_origin = this->_mesh.point(side->node(n)) - origin;

		    if(x_sym)
		      if( std::abs(dist_from_origin(0)) > 1.e-3 )
			on_x_sym=false;

		    if(y_sym)
		      if( std::abs(dist_from_origin(1)) > 1.e-3 )
			on_y_sym=false;

		    if(z_sym)
		      if( std::abs(dist_from_origin(2)) > 1.e-3 )
			on_z_sym=false;
	      
		    // 	      if(x_sym)
		    // 		if( std::abs(dist_from_origin(0)) > 1.e-6 )
		    // 		  on_x_sym=false;

		    // 	      if(y_sym)
		    // 		if( std::abs(dist_from_origin(1)) > 1.e-6 )
		    // 		  on_y_sym=false;

		    // 	      if(z_sym)
		    // 		if( std::abs(dist_from_origin(2)) > 1.e-6 )
		    // 		  on_z_sym=false;

		    //find the node most distant from origin
	      
		    Real r = dist_from_origin.size();
		    if (r > max_r)
		      {
			max_r = r;
			max_r_node=side->node(n);
		      }

		  }

		sym_side = (x_sym && on_x_sym) || (y_sym && on_y_sym) || (z_sym && on_z_sym);
				
		if (!sym_side)
		  faces.insert( std::make_pair(elem->id(), s) );					
					    
	      } // neighbor(s) == NULL
	  } // sides
      } // elems
  }



  
			

  //  If a boundary side has one node on the outer boundary,
  //  all points of this side are on the outer boundary.
  //  Start with the node most distant from origin, which has
  //  to be on the outer boundary, then recursively find all
  //  sides and nodes connected to it. Found sides are moved
  //  from faces to ofaces, nodes are collected in onodes.  
  //  Here, the search is done iteratively, because, depending on
  //  the mesh, a very high level of recursion might be necessary.
  onodes.insert(max_r_node);	


  {
    std::set< std::pair<unsigned int,unsigned int> >::iterator face_it = faces.begin();
    unsigned int facesfound=0;
    do {
			
      std::pair<unsigned int, unsigned int> p;
      p = *face_it;

      // This has to be a full-ordered side element,
      // since we need the correct n_nodes,
      AutoPtr<Elem> side(this->_mesh.elem(p.first)->build_side(p.second));
		
      bool found=false;
      for(unsigned int sn=0; sn<side->n_nodes(); sn++)
	if(onodes.count(side->node(sn)))
	  {
	    found=true;
	    break;
	  }
    	
    	
      // If a new oface is found, include its nodes in onodes
      if(found)		
	{
	  for(unsigned int sn=0; sn<side->n_nodes(); sn++)
	    onodes.insert(side->node(sn));
    				
	  ofaces.insert(p);
	  face_it++;			// iteration is done here
	  faces.erase(p);
    		
	  facesfound++;
	}
    
      else
	face_it++;			// iteration is done here


      // If at least one new oface was found in this cycle, 
      // do another search cycle.
      if(facesfound>0 && face_it == faces.end())
	{	
	  facesfound = 0;
	  face_it    = faces.begin();
	}
    
    } while(face_it != faces.end());
  }
  
	
#ifdef DEBUG
  if (be_verbose)
    std::cout << "  found " 
	      << faces.size() 
	      << " inner and " 
	      << ofaces.size() 
	      << " outer boundary faces" 
	      << std::endl;
#endif	
	
  // When the user provided a non-null pointer to
  // inner_faces, that implies he wants to have
  // this std::set.  For now, simply copy the data.
  if (inner_faces != NULL)
    *inner_faces = faces;

  // free memory, clear our local variable, no need
  // for it any more.
  faces.clear();


  // outer_nodes maps onodes to their duplicates
  std::map<unsigned int, Node *> outer_nodes;

  // We may need to pick our own object ids in parallel
  unsigned int old_max_node_id = _mesh.max_node_id();
  unsigned int old_max_elem_id = _mesh.max_elem_id();

  // for each boundary node, add an outer_node with 
  // double distance from origin.
  std::set<unsigned int>::iterator on_it = onodes.begin();
  for( ; on_it != onodes.end(); ++on_it)
    {
      Point p = (Point(this->_mesh.point(*on_it)) * 2.) - origin;
      if (_mesh.is_serial())
        {
          // Add with a default id in serial
          outer_nodes[*on_it]=this->_mesh.add_point(p);
        }
      else
        {
          // Pick a unique id in parallel
          Node &bnode = _mesh.node(*on_it);
          unsigned int new_id = bnode.id() + old_max_node_id;
          outer_nodes[*on_it] = 
            this->_mesh.add_point(p, new_id,
                                  bnode.processor_id());
        }
    }


#ifdef DEBUG
  // for verbose, remember n_elem
  unsigned int n_conventional_elem = this->_mesh.n_elem();
#endif


  // build Elems based on boundary side type
  std::set< std::pair<unsigned int,unsigned int> >::iterator face_it = ofaces.begin();
  for( ; face_it != ofaces.end(); ++face_it)
    {
      // Shortcut to the pair being iterated over
      std::pair<unsigned int,unsigned int> p = *face_it;
	
      // build a full-ordered side element to get the base nodes
      AutoPtr<Elem> side(this->_mesh.elem(p.first)->build_side(p.second)); 

      // create cell depending on side type, assign nodes,
      // use braces to force scope.
      bool is_higher_order_elem = false;
      {
        Elem* el;
	switch(side->type())
	{
	  // 3D infinite elements
	  // TRIs					
	  case TRI3:	
	    el=new InfPrism6;
	    break;
					 		
	  case TRI6: 
	    el=new InfPrism12; 
	    is_higher_order_elem = true;
	    break;
							
	  // QUADs					
	  case QUAD4: 
	    el=new InfHex8;
	    break;
							
	  case QUAD8: 
	    el=new InfHex16;
	    is_higher_order_elem = true;
	    break;
							
	  case QUAD9: 
	    el=new InfHex18;

	    // the method of assigning nodes (which follows below)
	    // omits in the case of QUAD9 the bubble node; therefore
	    // we assign these first by hand here.
	    el->set_node(16) = side->get_node(8);
	    el->set_node(17) = outer_nodes[side->node(8)];
	    is_higher_order_elem=true;
	    break;

	  // 2D infinite elements		 
	  case EDGE2:
	    el=new InfQuad4;
	    break;

	  case EDGE3:
	    el=new InfQuad6;
	    el->set_node(4) = side->get_node(2);
	    break;

	  // 1D infinite elements not supported
	  default: 
	    std::cout << "InfElemBuilder::build_inf_elem(Point, bool, bool, bool, bool): "
		      << "invalid face element "
		      << std::endl;
	    continue;
	}

        // In parallel, assign unique ids to the new element
        if (!_mesh.is_serial())
          {
            Elem *belem = _mesh.elem(p.first);
            el->processor_id() = belem->processor_id();
            // We'd better not have elements with more than 6 sides
            el->set_id (belem->id() * 6 + p.second + old_max_elem_id);
          }

	// assign vertices to the new infinite element
	const unsigned int n_base_vertices = side->n_vertices();
	for(unsigned int i=0; i<n_base_vertices; i++)
	  {
	    el->set_node(i                ) = side->get_node(i);
	    el->set_node(i+n_base_vertices) = outer_nodes[side->node(i)];
	  }


	// when this is a higher order element, 
	// assign also the nodes in between
	if (is_higher_order_elem)
          {
	    // n_safe_base_nodes is the number of nodes in \p side
	    // that may be safely assigned using below for loop.
	    // Actually, n_safe_base_nodes is _identical_ with el->n_vertices(),
	    // since for QUAD9, the 9th node was already assigned above
	    const unsigned int n_safe_base_nodes   = el->n_vertices();

	    for(unsigned int i=n_base_vertices; i<n_safe_base_nodes; i++)
	      {
	        el->set_node(i+n_base_vertices)   = side->get_node(i);
	        el->set_node(i+n_safe_base_nodes) = outer_nodes[side->node(i)];
	      }
	  }
			

	// add infinite element to mesh
	this->_mesh.add_elem(el);
      } // el goes out of scope
    } // for


#ifdef DEBUG
  ParallelMesh *pmesh = dynamic_cast<ParallelMesh *>(&this->_mesh);
  if (pmesh)
    pmesh->libmesh_assert_valid_parallel_ids();

  if (be_verbose)
    std::cout << "  added "
	      << this->_mesh.n_elem() - n_conventional_elem
	      << " infinite elements and "
	      << onodes.size() 
	      << " nodes to the mesh"
	      << std::endl
	      << std::endl;
#endif

}





#endif // ENABLE_INFINITE_ELEMENTS