dirichletbc.cpp

利用C

// Copyright (C) 2007-2008 Anders Logg and Garth N. Wells.
// Licensed under the GNU LGPL Version 2.1.
//
// Modified by Kristian Oelgaard, 2007
//
// First added:  2007-04-10
// Last changed: 2008-05-22

#include <dolfin/common/constants.h>
#include <dolfin/log/log.h>
#include <dolfin/mesh/Mesh.h>
#include <dolfin/mesh/MeshData.h>
#include <dolfin/mesh/Vertex.h>
#include <dolfin/mesh/Cell.h>
#include <dolfin/mesh/Facet.h>
#include <dolfin/mesh/Point.h>
#include <dolfin/mesh/SubDomain.h>
#include <dolfin/la/GenericMatrix.h>
#include <dolfin/la/GenericVector.h>
#include "Form.h"
#include "UFCMesh.h"
#include "UFCCell.h"
#include "SubSystem.h"
#include "DirichletBC.h"

using namespace dolfin;

//-----------------------------------------------------------------------------
DirichletBC::DirichletBC(Function& g,
                         Mesh& mesh,
                         SubDomain& sub_domain,
                         BCMethod method)
  : BoundaryCondition(), g(g), _mesh(mesh),
    method(method), user_sub_domain(&sub_domain)
{
  initFromSubDomain(sub_domain);
}
//-----------------------------------------------------------------------------
DirichletBC::DirichletBC(Function& g,
                         MeshFunction<uint>& sub_domains,
                         uint sub_domain,
                         BCMethod method)
  : BoundaryCondition(), g(g), _mesh(sub_domains.mesh()),
    method(method), user_sub_domain(0)
{
  initFromMeshFunction(sub_domains, sub_domain);
}
//-----------------------------------------------------------------------------
DirichletBC::DirichletBC(Function& g,
                         Mesh& mesh,
                         uint sub_domain,
                         BCMethod method)
  : BoundaryCondition(), g(g), _mesh(mesh),
    method(method), user_sub_domain(0)
{
  initFromMesh(sub_domain);
}
//-----------------------------------------------------------------------------
DirichletBC::DirichletBC(Function& g,
                         Mesh& mesh,
                         SubDomain& sub_domain,
                         const SubSystem& sub_system,
                         BCMethod method)
  : BoundaryCondition(), g(g), _mesh(mesh),
    method(method), user_sub_domain(&sub_domain),
    sub_system(sub_system)
{
  initFromSubDomain(sub_domain);
}
//-----------------------------------------------------------------------------
DirichletBC::DirichletBC(Function& g,
                         MeshFunction<uint>& sub_domains,
                         uint sub_domain,
                         const SubSystem& sub_system,
                         BCMethod method)
  : BoundaryCondition(), g(g), _mesh(sub_domains.mesh()),
    method(method), user_sub_domain(0),
    sub_system(sub_system)
{
  initFromMeshFunction(sub_domains, sub_domain);
}
//-----------------------------------------------------------------------------
DirichletBC::DirichletBC(Function& g,
                         Mesh& mesh,
                         uint sub_domain,
                         const SubSystem& sub_system,
                         BCMethod method)
  : BoundaryCondition(), g(g), _mesh(mesh),
    method(method), user_sub_domain(0),
    sub_system(sub_system)
{
  initFromMesh(sub_domain);
}
//-----------------------------------------------------------------------------
DirichletBC::~DirichletBC()
{
  // Do nothing
}
//-----------------------------------------------------------------------------
void DirichletBC::apply(GenericMatrix& A, GenericVector& b, const Form& form)
{
  apply(A, b, 0, form.dofMaps()[1], form.form());
}
//-----------------------------------------------------------------------------
void DirichletBC::apply(GenericMatrix& A, GenericVector& b, const DofMap& dof_map,
                        const ufc::form& form)
{
  apply(A, b, 0, dof_map, form);
}
//-----------------------------------------------------------------------------
void DirichletBC::apply(GenericMatrix& A, GenericVector& b,
                        const GenericVector& x, const Form& form)
{
  apply(A, b, &x, form.dofMaps()[1], form.form());
}
//-----------------------------------------------------------------------------
void DirichletBC::apply(GenericMatrix& A, GenericVector& b,
                        const GenericVector& x, const DofMap& dof_map, const ufc::form& form)
{
  apply(A, b, &x, dof_map, form);
}
//-----------------------------------------------------------------------------
void DirichletBC::apply(GenericMatrix& A, GenericVector& b,
                        const GenericVector* x, const DofMap& dof_map, const ufc::form& form)
{
  // Simple check
  const uint N = dof_map.global_dimension();
  if (N != A.size(0) /*  || N != A.size(1) alow for rectangular matrices */)
    error("Incorrect dimension of matrix for application of boundary conditions. Did you assemble it on a different mesh?");
  if (N != b.size())
    error("Incorrect dimension of matrix for application of boundary conditions. Did you assemble it on a different mesh?");

  // A map to hold the mapping from boundary dofs to boundary values
  std::map<uint, real> boundary_values;

  // Create local data for application of boundary conditions
  BoundaryCondition::LocalData data(form, _mesh, dof_map, sub_system);

  // Compute dofs and values
  computeBC(boundary_values, data);

  // Copy boundary value data to arrays
  uint* dofs = new uint[boundary_values.size()];
  real* values = new real[boundary_values.size()];
  std::map<uint, real>::const_iterator boundary_value;
  uint i = 0;
  for (boundary_value = boundary_values.begin(); boundary_value != boundary_values.end(); ++boundary_value)
  {
    dofs[i]     = boundary_value->first;
    values[i++] = boundary_value->second;
  }
  
  // Modify boundary values for nonlinear problems
  if (x)
  {
    real* x_values = new real[boundary_values.size()];
    x->get(x_values, boundary_values.size(), dofs);
    for (uint i = 0; i < boundary_values.size(); i++)
      values[i] -= x_values[i];
  }
  
  message("Applying boundary conditions to linear system.");
  
  // Modify RHS vector (b[i] = value)
  b.set(values, boundary_values.size(), dofs);
  
  // Modify linear system (A_ii = 1)
  A.ident(boundary_values.size(), dofs);
  
  // Clear temporary arrays
  delete [] dofs;
  delete [] values;
  
  // Finalise changes to A
  A.apply();
  
  // Finalise changes to b
  b.apply();
}
//-----------------------------------------------------------------------------
void DirichletBC::zero(GenericMatrix& A, const Form& form)
{
  zero(A, form.dofMaps()[1], form.form());
}
//-----------------------------------------------------------------------------
void DirichletBC::zero(GenericMatrix& A, const DofMap& dof_map, const ufc::form& form)
{
  // Simple check
  const uint N = dof_map.global_dimension();
  if (N != A.size(0))
    error("Incorrect dimension of matrix for application of boundary conditions. Did you assemble it on a different mesh?");

  // A map to hold the mapping from boundary dofs to boundary values
  std::map<uint, real> boundary_values;

  // Create local data for application of boundary conditions
  BoundaryCondition::LocalData data(form, _mesh, dof_map, sub_system);

  // Compute dofs and values
  computeBC(boundary_values, data);

  // Copy boundary value data to arrays
  uint* dofs = new uint[boundary_values.size()];
  std::map<uint, real>::const_iterator boundary_value;
  uint i = 0;
  for (boundary_value = boundary_values.begin(); boundary_value != boundary_values.end(); ++boundary_value)
    dofs[i++] = boundary_value->first;

  // Modify linear system (A_ii = 1)
  A.zero(boundary_values.size(), dofs);

  // Finalise changes to A
  A.apply();

  // Clear temporary arrays
  delete [] dofs;
}
//-----------------------------------------------------------------------------
Mesh& DirichletBC::mesh()
{
  return _mesh;
}
//-----------------------------------------------------------------------------
void DirichletBC::initFromSubDomain(SubDomain& sub_domain)
{
  dolfin_assert(facets.size() == 0);

  // FIXME: This can be made more efficient, we should be able to
  // FIXME: extract the facets without first creating a MeshFunction on
  // FIXME: the entire mesh and then extracting the subset. This is done
  // FIXME: mainly for convenience (we may reuse mark() in SubDomain).

  // Make sure the mesh has been ordered
  _mesh.order();

  // Create mesh function for sub domain markers on facets
  const uint dim = _mesh.topology().dim();
  _mesh.init(dim - 1);
  MeshFunction<uint> sub_domains(_mesh, dim - 1);

  // Mark everything as sub domain 1
  sub_domains = 1;
  
  // Mark the sub domain as sub domain 0
  sub_domain.mark(sub_domains, 0);

  // Initialize from mesh function
  initFromMeshFunction(sub_domains, 0);
}
//-----------------------------------------------------------------------------
void DirichletBC::initFromMeshFunction(MeshFunction<uint>& sub_domains,
                                       uint sub_domain)
{
  dolfin_assert(facets.size() == 0);

  // Make sure we have the facet - cell connectivity
  const uint dim = _mesh.topology().dim();
  _mesh.init(dim - 1, dim);

  // Make sure the mesh has been ordered
  _mesh.order();

  // Build set of boundary facets
  for (FacetIterator facet(_mesh); !facet.end(); ++facet)
  {
    // Skip facets not on this boundary
    if (sub_domains(*facet) != sub_domain)
      continue;

    // Get cell to which facet belongs (there may be two, but pick first)
    Cell cell(_mesh, facet->entities(dim)[0]);