dirichletbc.cpp

利用C

    // Get local index of facet with respect to the cell
    const uint facet_number = cell.index(*facet);

    // Copy data
    facets.push_back(std::pair<uint, uint>(cell.index(), facet_number));
  }
}
//-----------------------------------------------------------------------------
void DirichletBC::initFromMesh(uint sub_domain)
{
  dolfin_assert(facets.size() == 0);

  // For this to work, the mesh *needs* to be ordered according to
  // the UFC ordering before it gets here. So reordering the mesh
  // here will either have no effect (if the mesh is already ordered
  // or it won't do anything good (since the markers are wrong anyway).
  // In conclusion: we don't need to order the mesh here.
  
  cout << "Creating sub domain markers for boundary condition." << endl;

  // Get data
  Array<uint>* facet_cells   = _mesh.data().array("boundary facet cells");
  Array<uint>* facet_numbers = _mesh.data().array("boundary facet numbers");
  Array<uint>* indicators    = _mesh.data().array("boundary indicators");

  // Check data
  if (!facet_cells)
    error("Mesh data \"boundary facet cells\" not available.");
  if (!facet_numbers)
    error("Mesh data \"boundary facet numbers\" not available.");
  if (!indicators)
    error("Mesh data \"boundary indicators\" not available.");

  // Get size
  const uint size = facet_cells->size();
  dolfin_assert(size == facet_numbers->size());
  dolfin_assert(size == indicators->size());

  // Build set of boundary facets
  for (uint i = 0; i < size; i++)
  {
    // Skip facets not on this boundary
    if ((*indicators)[i] != sub_domain)
      continue;

    // Copy data
    facets.push_back(std::pair<uint, uint>((*facet_cells)[i], (*facet_numbers)[i]));
  }
}
//-----------------------------------------------------------------------------
void DirichletBC::computeBC(std::map<uint, real>& boundary_values,
                            BoundaryCondition::LocalData& data)
{
  // Choose strategy
  switch (method)
  {
  case topological:
    computeBCTopological(boundary_values, data);
    break;
  case geometric:
    computeBCGeometric(boundary_values, data);
    break;
  case pointwise:
    computeBCPointwise(boundary_values, data);
    break;
  default:
    error("Unknown method for application of boundary conditions.");
  }
}
//-----------------------------------------------------------------------------
void DirichletBC::computeBCTopological(std::map<uint, real>& boundary_values,
                                       BoundaryCondition::LocalData& data)
{
  // Special case
  if (facets.size() == 0)
  {
    warning("Found no facets matching domain for boundary condition.");
    return;
  }

  // Iterate over facets
  Progress p("Computing Dirichlet boundary values, topological search", facets.size());
  for (uint f = 0; f < facets.size(); f++)
  {
    // Get cell number and local facet number
    uint cell_number = facets[f].first;
    uint facet_number = facets[f].second;

    // Create cell
    Cell cell(_mesh, cell_number);
    UFCCell ufc_cell(cell);

    // Interpolate function on cell
    g.interpolate(data.w, ufc_cell, *data.finite_element, cell, facet_number);
    
    // Tabulate dofs on cell
    data.dof_map->tabulate_dofs(data.cell_dofs, ufc_cell);
    
    // Tabulate which dofs are on the facet
    data.dof_map->tabulate_facet_dofs(data.facet_dofs, facet_number);
    
    // Debugging print:
    /* 
       cout << endl << "Handling BC's for:" << endl;
       cout << "Cell:  " << facet.entities(facet.dim() + 1)[0] << endl;
       cout << "Facet: " << local_facet << endl;
    */
    
    // Pick values for facet
    for (uint i = 0; i < data.dof_map->num_facet_dofs(); i++)
    {
      const uint dof = data.offset + data.cell_dofs[data.facet_dofs[i]];
      const real value = data.w[data.facet_dofs[i]];
      boundary_values[dof] = value;
      //cout << "Setting BC value: i = " << i << ", dof = " << dof << ", value = " << value << endl;
    }

    p++;
  }
}
//-----------------------------------------------------------------------------
void DirichletBC::computeBCGeometric(std::map<uint, real>& boundary_values,
                                     BoundaryCondition::LocalData& data)
{
  // Special case
  if (facets.size() == 0)
  {
    warning("Found no facets matching domain for boundary condition.");
    return;
  }

  // Initialize facets, needed for geometric search
  message("Computing facets, needed for geometric application of boundary conditions.");
  _mesh.init(_mesh.topology().dim() - 1);

  // Iterate over facets
  Progress p("Computing Dirichlet boundary values, geometric search", facets.size());
  for (uint f = 0; f < facets.size(); f++)
  {
    // Get cell number and local facet number
    uint cell_number = facets[f].first;
    uint facet_number = facets[f].second;

    // Create facet
    Cell cell(_mesh, cell_number);
    Facet facet(_mesh, cell.entities(_mesh.topology().dim() - 1)[facet_number]);

    // Loop the vertices associated with the facet
    for (VertexIterator vertex(facet); !vertex.end(); ++vertex)
    {
      // Loop the cells associated with the vertex
      for (CellIterator c(*vertex); !c.end(); ++c)
      {
        UFCCell ufc_cell(*c);
        
        // Interpolate function on cell
        g.interpolate(data.w, ufc_cell, *data.finite_element, *c);
        
        // Tabulate dofs on cell, and their coordinates
        data.dof_map->tabulate_dofs(data.cell_dofs, ufc_cell);
        data.dof_map->tabulate_coordinates(data.coordinates, ufc_cell);
        
        // Loop all dofs on cell
        for (uint i = 0; i < data.dof_map->local_dimension(); ++i)
        {
          // Check if the coordinates are on current facet and thus on boundary
          if (!onFacet(data.coordinates[i], facet))
            continue;
          
          // Set boundary value
          const uint dof = data.offset + data.cell_dofs[i];
          const real value = data.w[i];
          boundary_values[dof] = value;
        }
      }
    }
  }
}
//-----------------------------------------------------------------------------
void DirichletBC::computeBCPointwise(std::map<uint, real>& boundary_values,
                                     BoundaryCondition::LocalData& data)
{
  dolfin_assert(user_sub_domain);

  // Iterate over cells
  Progress p("Computing Dirichlet boundary values, pointwise search", facets.size());
  for (CellIterator cell(_mesh); !cell.end(); ++cell)
  {
    // Interpolate function on cell
    UFCCell ufc_cell(*cell);
    g.interpolate(data.w, ufc_cell, *data.finite_element, *cell);
    
    // Tabulate dofs on cell, and their coordinates
    data.dof_map->tabulate_dofs(data.cell_dofs, ufc_cell);
    data.dof_map->tabulate_coordinates(data.coordinates, ufc_cell);
    
    // Loop all dofs on cell
    for (uint i = 0; i < data.dof_map->local_dimension(); ++i)
    {
      // Check if the coordinates are part of the sub domain
      if ( !user_sub_domain->inside(data.coordinates[i], false) )
        continue;
      
      // Set boundary value
      const uint dof = data.offset + data.cell_dofs[i];
      const real value = data.w[i];
      boundary_values[dof] = value;
    }

    p++;
  }
}
//-----------------------------------------------------------------------------
bool DirichletBC::onFacet(real* coordinates, Facet& facet)
{
  // Check if the coordinates are on the same line as the line segment
  if ( facet.dim() == 1 )
  {
    // Create points
    Point p(coordinates[0], coordinates[1]);
    Point v0 = Vertex(facet.mesh(), facet.entities(0)[0]).point();
    Point v1 = Vertex(facet.mesh(), facet.entities(0)[1]).point();

    // Create vectors
    Point v01 = v1 - v0;
    Point vp0 = v0 - p;
    Point vp1 = v1 - p;

    // Check if the length of the sum of the two line segments vp0 and vp1 is
    // equal to the total length of the facet
    if ( std::abs(v01.norm() - vp0.norm() - vp1.norm()) < DOLFIN_EPS )
      return true;
    else
      return false;
  }
  // Check if the coordinates are in the same plane as the triangular facet
  else if ( facet.dim() == 2 )
  {
    // Create points
    Point p(coordinates[0], coordinates[1], coordinates[2]);
    Point v0 = Vertex(facet.mesh(), facet.entities(0)[0]).point();
    Point v1 = Vertex(facet.mesh(), facet.entities(0)[1]).point();
    Point v2 = Vertex(facet.mesh(), facet.entities(0)[2]).point();

    // Create vectors
    Point v01 = v1 - v0;
    Point v02 = v2 - v0;
    Point vp0 = v0 - p;
    Point vp1 = v1 - p;
    Point vp2 = v2 - p;

    // Check if the sum of the area of the sub triangles is equal to the total
    // area of the facet
    if ( std::abs(v01.cross(v02).norm() - vp0.cross(vp1).norm() - vp1.cross(vp2).norm()
        - vp2.cross(vp0).norm()) < DOLFIN_EPS )
      return true;
    else
      return false;
  }

  error("Unable to determine if given point is on facet (not implemented for given facet dimension).");

  return false;
}
//-----------------------------------------------------------------------------
void DirichletBC::setSubSystem(SubSystem sub_system)
{
  this->sub_system = sub_system;
}
//-----------------------------------------------------------------------------