fe_lagrange.c

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

            return 4;

          default:
            {
#ifdef DEBUG
              std::cerr << "ERROR: Bad ElemType = " << t
                << " for THIRD order approximation!" 
                << std::endl;
#endif
              libmesh_error();	
            }
        }
      }
      
    default:
      libmesh_error();
    }
  
  libmesh_error();  
  return 0;
}



template <unsigned int Dim, FEFamily T>
unsigned int FE<Dim,T>::n_dofs_at_node(const ElemType t,
				       const Order o,
				       const unsigned int n)
{
  switch (o)
    {
      
      // linear Lagrange shape functions
    case FIRST:
      {
	switch (t)
	  {
	  case NODEELEM:
	    return 1;

	  case EDGE2:
	  case EDGE3:
          case EDGE4:
	    {
	      switch (n)
		{
		case 0:
		case 1:
		  return 1;
		  
		default:
		  return 0;
		}
	    }

	  case TRI3:
	  case TRI6:
	    {
	      switch (n)
		{
		case 0:
		case 1:
		case 2:
		  return 1;
		  
		default:
		  return 0;
		}
	    }

	  case QUAD4:
	  case QUAD8:
	  case QUAD9:
	    {
	      switch (n)
		{
		case 0:
		case 1:
		case 2:
		case 3:
		  return 1;
		  
		default:
		  return 0;
		}
	    }


	  case TET4:
	  case TET10:
	    {
	      switch (n)
		{
		case 0:
		case 1:
		case 2:
		case 3:
		  return 1;
		  
		default:
		  return 0;
		}
	    }
		    
	  case HEX8:
	  case HEX20:
	  case HEX27:
	    {
	      switch (n)
		{
		case 0:
		case 1:
		case 2:
		case 3:
		case 4:
		case 5:
		case 6:
		case 7:
		  return 1;
		  
		default:
		  return 0;
		}
	    }

	  case PRISM6:
	  case PRISM15:
	  case PRISM18:
	    {
	      switch (n)
		{
		case 0:
		case 1:
		case 2:
		case 3:
		case 4:
		case 5:
		  return 1;
		  
		default:
		  return 0;
		}
	    }

	  case PYRAMID5:
	    {
	      switch (n)
		{
		case 0:
		case 1:
		case 2:
		case 3:
		case 4:
		  return 1;
		  
		default:
		  return 0;
		}
	    }

	  default:
	    {
#ifdef DEBUG
	      std::cerr << "ERROR: Bad ElemType = " << t
			<< " for FIRST order approximation!" 
			<< std::endl;
#endif
	      libmesh_error();	    
	    }
	  }
      }

      // quadratic Lagrange shape functions
    case SECOND:
      {
	switch (t)
	  {
	  // quadratic lagrange has one dof at each node
	  case NODEELEM:
	  case EDGE3:
	  case TRI6:
	  case QUAD8:
	  case QUAD9:
	  case TET10:
	  case HEX20:
	  case HEX27:
	  case PRISM15:
	  case PRISM18:
	    return 1;

	  default:
	    {
#ifdef DEBUG
	      std::cerr << "ERROR: Bad ElemType = " << t
			<< " for SECOND order approximation!" 
			<< std::endl;
#endif
	      libmesh_error();	    
	    }
	  }
      }

    case THIRD:
      {
        switch (t)
          {
	    case NODEELEM:
            case EDGE4:
              return 1;

	  default:
	    {
#ifdef DEBUG
	      std::cerr << "ERROR: Bad ElemType = " << t
			<< " for THIRD order approximation!" 
			<< std::endl;
#endif
	      libmesh_error();	    
	    }
          }
      }


      
    default:
      libmesh_error();
    }
  
  libmesh_error();  
  return 0;
}



template <unsigned int Dim, FEFamily T>
unsigned int FE<Dim,T>::n_dofs_per_elem(const ElemType,
					const Order)
{
  // Lagrange elements have no dofs per element
  // (just at the nodes)
  
  return 0;
}



template <unsigned int Dim, FEFamily T>
FEContinuity FE<Dim,T>::get_continuity() const
{
  return C_ZERO;
}



template <unsigned int Dim, FEFamily T>
bool FE<Dim,T>::is_hierarchic() const
{
  return false;
}



#ifdef ENABLE_AMR
template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::compute_constraints (DofConstraints &constraints,
				     DofMap &dof_map,
				     const unsigned int variable_number,
				     const Elem* elem)
{
  // Only constrain elements in 2,3D.
  if (Dim == 1)
    return;

  libmesh_assert (elem != NULL);

  // Only constrain active and ancestor elements
  if (elem->subactive())
    return;

  FEType fe_type = dof_map.variable_type(variable_number);
  fe_type.order = static_cast<Order>(fe_type.order + elem->p_level());

  std::vector<unsigned int> my_dof_indices, parent_dof_indices;

  // Look at the element faces.  Check to see if we need to 
  // build constraints.
  for (unsigned int s=0; s<elem->n_sides(); s++)
    if (elem->neighbor(s) != NULL)
      if (elem->neighbor(s)->level() < elem->level()) // constrain dofs shared between
	{                                                     // this element and ones coarser
	                                                      // than this element.
	  // Get pointers to the elements of interest and its parent.
	  const Elem* parent = elem->parent();
	  
	  // This can't happen...  Only level-0 elements have NULL
	  // parents, and no level-0 elements can be at a higher
	  // level than their neighbors!
	  libmesh_assert (parent != NULL);
	  
	  const AutoPtr<Elem> my_side     (elem->build_side(s));
	  const AutoPtr<Elem> parent_side (parent->build_side(s));

	  // This function gets called element-by-element, so there
	  // will be a lot of memory allocation going on.  We can 
	  // at least minimize this for the case of the dof indices
	  // by efficiently preallocating the requisite storage.
	  my_dof_indices.reserve (my_side->n_nodes());
	  parent_dof_indices.reserve (parent_side->n_nodes());

	  dof_map.dof_indices (my_side.get(), my_dof_indices,
			       variable_number);
	  dof_map.dof_indices (parent_side.get(), parent_dof_indices,
			       variable_number);
  
	  for (unsigned int my_dof=0;
	       my_dof<FEInterface::n_dofs(Dim-1, fe_type, my_side->type());
	       my_dof++)
	    {
	      libmesh_assert (my_dof < my_side->n_nodes());
	      
	      // My global dof index.
	      const unsigned int my_dof_g = my_dof_indices[my_dof];
	      
	      // The support point of the DOF
	      const Point& support_point = my_side->point(my_dof);
	      
	      // Figure out where my node lies on their reference element.
	      const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
								  parent_side.get(),
								  support_point);
	      
	      // Compute the parent's side shape function values.
	      for (unsigned int their_dof=0;
		   their_dof<FEInterface::n_dofs(Dim-1, fe_type, parent_side->type());
		   their_dof++)
		{
		  libmesh_assert (their_dof < parent_side->n_nodes());
		  
	          // Their global dof index.
	          const unsigned int their_dof_g =
				  parent_dof_indices[their_dof];
		  
		  const Real their_dof_value = FEInterface::shape(Dim-1,
								  fe_type,
								  parent_side->type(),
								  their_dof,
								  mapped_point);
		  
		  // Only add non-zero and non-identity values
		  // for Lagrange basis functions.
		  if ((std::abs(their_dof_value) > 1.e-5) &&
		      (std::abs(their_dof_value) < .999)) 
		    {
		      // since we may be running this method concurretly 
		      // on multiple threads we need to acquire a lock 
		      // before modifying the shared constraint_row object.		      
		      Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);

		      // A reference to the constraint row.
		      DofConstraintRow& constraint_row = constraints[my_dof_g];
		      
		      constraint_row.insert(std::make_pair (their_dof_g,
							    their_dof_value));
		    }
#ifdef DEBUG
		  // Protect for the case u_i = 0.999 u_j,
		  // in which case i better equal j.
		  else if (their_dof_value >= .999)
		    libmesh_assert (my_dof_g == their_dof_g);
#endif
		}		      
	    }
	}
}
#endif // #ifdef ENABLE_AMR



template <unsigned int Dim, FEFamily T>
bool FE<Dim,T>::shapes_need_reinit() const
{
  return false;
}




//--------------------------------------------------------------
// Explicit instantiation of member functions
INSTANTIATE_MBRF(1,LAGRANGE);
INSTANTIATE_MBRF(2,LAGRANGE);
INSTANTIATE_MBRF(3,LAGRANGE);

#ifdef ENABLE_AMR
template void FE<2,LAGRANGE>::compute_constraints(DofConstraints&, DofMap&, 
						  const unsigned int,
						  const Elem*);
template void FE<3,LAGRANGE>::compute_constraints(DofConstraints&, DofMap&, 
						  const unsigned int,
						  const Elem*);
#endif // #ifdef ENABLE_AMR