[Getfem-users] Dirichlet Boundary Conditions

[John Burnell]

Hi

I have been trying to setup a few problems in getfem, but I am having 
difficulties with dirichlet boundary conditions. I have modified the 
Stokes example from the tests directory to show the problem.

Basically I want to have different  Dirichlet conditions on different 
parts of the boundary (a velocity prescribed on 1 part and noflow on 
others). So I setup 2 different regions and applied the boundary 
conditions. But I end up with a singular matrix in SuperLU.

In the test I attach, I started with the Stokes example and divided the 
Dirichlet boundary into 2 portions. I then set velocities equal to 0 on 
each portion (so the mathematical problem is the same as the Stokes 
example). When I run the example I get the error "SuperLU solve failed: 
info = 2216". My new code is marked with !!!NEW.

Does anyone understand what is going on?

I can resolve the problem by setting one of the  Dirichlet bricks to 
"Penalization", but this isn't entirely satisfactory.

Thanks
John Burnell

--
==========================================
John G. Burnell
Research Scientist
Applied Mathematics Team
Industrial Research
P O Box 31-310, Lower Hutt
New Zealand.                       
Telephone: 64-4-9313-352       
Fax:       64-4-9313-003        
EMail:  address@hidden


This message has been scanned for viruses and dangerous content by MailScanner, 
and is believed to be clean.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parameters for program stokes                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% pde parameters :                                                %%%%%
LX = 0.2e-3;            % size in X.
LY = 0.4e-3;            % size in Y.
LZ = 1.0e-3;            % size in Z.
NU = 1.0;               % Viscosity coefficient.

%%%%%   discretisation parameters  :                                  %%%%%
MESH_TYPE = 'GT_PK(3,1)';         % linear tetrahedra
%MESH_TYPE = 'GT_LINEAR_QK(2,1)'; % linear rectangles
%MESH_TYPE = 'GT_PRISM(3,1)';     % 3D prisms
NX = 3;                           % space step.
MESH_NOISED = 0; % Set to one if you want to "shake" the mesh

FEM_TYPE = 'FEM_PK_WITH_CUBIC_BUBBLE(3, 2)';
%FEM_TYPE = 'FEM_PK(3,3)';  % PK method
%FEM_TYPE = 'FEM_QK(2,1)';  % QK method
%FEM_TYPE = 'FEM_PRODUCT(FEM_PK(2,1),FEM_PK(1,1))'; % tensorial product of FEM 
for prisms
%FEM_TYPE = 'FEM_PK_HIERARCHICAL(2,2)'; % Hierarchical PK on simplexes
%FEM_TYPE = 'FEM_PK_HIERARCHICAL_COMPOSITE(2,1,2)'; % Hierarchical PK with s 
divisions

FEM_TYPE_P = 'FEM_PK(3,1)';  % P1 for triangles
% FEM_TYPE_P = 'FEM_PK_DISCONTINUOUS(3,1)';  % Discontinuous P1 for triangles

% DATA_FEM_TYPE must be defined if your main FEM is not Lagrangian
DATA_FEM_TYPE = 'FEM_PK(3,1)';

%INTEGRATION = 'IM_TRIANGLE(6)'; % quadrature rule for polynomials up
                               % to degree 6 on triangles
INTEGRATION = 'IM_TETRAHEDRON(6)'; % quadrature rule for polynomials up
                               % to degree 6 on tetra
%INTEGRATION = 'IM_EXACT_SIMPLEX(2)'; % exact integration on triangles
%INTEGRATION = 'IM_NC(2,6)';     % newton-cotes of degree 6 on triangles
%INTEGRATION = 'IM_NC_PARALLELEPIPED(2,6)'; % newton-cotes, degree 6,
                                          % quadrangles
%INTEGRATION = 'IM_NC_PRISM(3,12)'; % newton-cotes, degree 12, prims
%INTEGRATION = 'IM_GAUSS1D(10)'; % Gauss-Legendre integration on the
                               % segment of order 10
%INTEGRATION = 'IM_GAUSSLOBATTO1D(10)'; % Gauss-Lobatto-Legendre
                                      % integration on the segment
                                      % of order 10
%INTEGRATION = 'IM_GAUSS_PARALLELEPIPED(2,10)'; % Product of two
                                              % IM_GAUSS1D(10) (for
                                              % quadrangles)
%INTEGRATION = 'IM_STRUCTURED_COMPOSITE(IM_GAUSS1D(5), 3)';
%INTEGRATION = 'IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(7), 3)';

RESIDUAL = 1E-9;        % residu for conjugate gradient.

%%%%%   saving parameters                                             %%%%%
ROOTFILENAME = 'stokes';     % Root of data files.
VTK_EXPORT = 1 % export solution to a .vtk file ?

// -*- c++ -*- (enables emacs c++ mode)

#include <getfem_assembling.h> /* import assembly methods (and norms comp.) */
#include <getfem_export.h>   /* export functions (save solution in a file)  */
#include <getfem_regular_meshes.h>
#include <getfem_model_solvers.h>
#include <gmm.h>

/* some Getfem++ types that we will be using */
using bgeot::base_small_vector; /* special class for small (dim<16) vectors */
using bgeot::base_node;  /* geometrical nodes(derived from base_small_vector)*/
using bgeot::scalar_type; /* = double */
using bgeot::size_type;   /* = unsigned long */
using bgeot::base_matrix; /* small dense matrix. */

/* definition of some matrix/vector types. These ones are built
 * using the predefined types in Gmm++
 */
typedef getfem::modeling_standard_sparse_vector sparse_vector;
typedef getfem::modeling_standard_sparse_matrix sparse_matrix;
typedef getfem::modeling_standard_plain_vector  plain_vector;

/*
 * structure for the stokes problem
 */
struct stokes_problem {

  // !!!NEW added LEFT_BOUNDARY_NUM
  enum { DIRICHLET_BOUNDARY_NUM = 0,
         NEUMANN_BOUNDARY_NUM,
         LEFT_BOUNDARY_NUM
  };
  getfem::mesh mesh;         /* the mesh */
  getfem::mesh_im  mim;      /* integration methods.                         */
  getfem::mesh_fem mf_u;     /* main mesh_fem, for the velocity              */
  getfem::mesh_fem mf_p;     /* mesh_fem for the pressure                    */
  getfem::mesh_fem mf_rhs;   /* mesh_fem for the right hand side (f(x),..)   */
  scalar_type nu;            /* Lamé coefficients.                           */

  scalar_type residual;        /* max residual for the iterative solvers        
 */

  std::string datafilename;
  ftool::md_param PARAM;

  bool solve(plain_vector &U);
  void init(void);
  stokes_problem(void) : mim(mesh), mf_u(mesh), mf_p(mesh),
                         mf_rhs(mesh) {}
};

/* Read parameters from the .param file, build the mesh, set finite element
 * and integration methods and selects the boundaries.
 */
void stokes_problem::init(void) {
  std::string MESH_TYPE = PARAM.string_value("MESH_TYPE","Mesh type ");
  std::string FEM_TYPE  = PARAM.string_value("FEM_TYPE","FEM name");
  std::string FEM_TYPE_P  = PARAM.string_value("FEM_TYPE_P","FEM name P");
  std::string INTEGRATION = PARAM.string_value("INTEGRATION",
                                               "Name of integration method");
  cout << "MESH_TYPE=" << MESH_TYPE << "\n";
  cout << "FEM_TYPE="  << FEM_TYPE << "\n";
  cout << "INTEGRATION=" << INTEGRATION << "\n";

  /* First step : build the mesh */
  bgeot::pgeometric_trans pgt = 
    bgeot::geometric_trans_descriptor(MESH_TYPE);
  size_type N = pgt->dim();
  std::vector<size_type> nsubdiv(N);
  std::fill(nsubdiv.begin(),nsubdiv.end(),
            PARAM.int_value("NX", "Nomber of space steps "));
  getfem::regular_unit_mesh(mesh, nsubdiv, pgt,
                            PARAM.int_value("MESH_NOISED") != 0);
  
  /* scale the unit mesh to [LX,LY,..] and incline it */
   bgeot::base_matrix M(N,N);
  for (size_type i=0; i < N; ++i) {
    static const char *t[] = {"LX","LY","LZ"};
    M(i,i) = (i<3) ? PARAM.real_value(t[i],t[i]) : 1.0;
  }
  mesh.transformation(M);

  datafilename = PARAM.string_value("ROOTFILENAME","Base name of data files.");
  residual = PARAM.real_value("RESIDUAL"); if (residual == 0.) residual = 1e-10;

  nu = PARAM.real_value("NU", "Viscosité");
  mf_u.set_qdim(N);

  /* set the finite element on the mf_u */
  getfem::pfem pf_u = 
    getfem::fem_descriptor(FEM_TYPE);
  getfem::pintegration_method ppi = 
    getfem::int_method_descriptor(INTEGRATION); 

  mim.set_integration_method(mesh.convex_index(), ppi);
  mf_u.set_finite_element(mesh.convex_index(), pf_u);
  mf_p.set_finite_element(mesh.convex_index(),
                          getfem::fem_descriptor(FEM_TYPE_P));

  /* set the finite element on mf_rhs (same as mf_u is DATA_FEM_TYPE is
     not used in the .param file */
  std::string data_fem_name = PARAM.string_value("DATA_FEM_TYPE");
  if (data_fem_name.size() == 0) {
    if (!pf_u->is_lagrange()) {
      DAL_THROW(dal::failure_error, "You are using a non-lagrange FEM "
                << data_fem_name << ". In that case you need to set "
                << "DATA_FEM_TYPE in the .param file");
    }
    mf_rhs.set_finite_element(mesh.convex_index(), pf_u);
  } else {
    mf_rhs.set_finite_element(mesh.convex_index(), 
                              getfem::fem_descriptor(data_fem_name));
  }
  
  /* set boundary conditions
   * (Neuman on the upper face, Dirichlet elsewhere) */
  cout << "Selecting Neumann and Dirichlet boundaries\n";
  getfem::mesh_region border_faces;
  getfem::outer_faces_of_mesh(mesh, border_faces);
  for (getfem::mr_visitor it(border_faces); !it.finished(); ++it) {
    assert(it.is_face());
    base_node un = mesh.normal_of_face_of_convex(it.cv(), it.f());
    un /= gmm::vect_norm2(un);
    if (gmm::abs(un[N-1] - 1.0) < 1.0E-7) { // new Neumann face
      mesh.region(NEUMANN_BOUNDARY_NUM).add(it.cv(), it.f());
    } else if (gmm::abs(un[0] - 1.0) < 1.0E-7) { // !!!NEW left boundary
      mesh.region(LEFT_BOUNDARY_NUM).add(it.cv(), it.f());
    } else {
      mesh.region(DIRICHLET_BOUNDARY_NUM).add(it.cv(), it.f());
    }
  }
}

/**************************************************************************/
/*  Model.                                                                */
/**************************************************************************/

base_small_vector sol_f(const base_small_vector &P) {
  base_small_vector res(P.size());
  res[P.size()-1] = -10.0;
  return res;
}


bool stokes_problem::solve(plain_vector &U) {
  size_type nb_dof_rhs = mf_rhs.nb_dof();
  size_type N = mesh.dim();

  // Linearized elasticity brick.
  getfem::mdbrick_isotropic_linearized_elasticity<>
    ELAS(mim, mf_u, 0.0, nu);

  // Pressure term
  getfem::mdbrick_linear_incomp<> INCOMP(ELAS, mf_p);
  
  // Defining the volumic source term.
  plain_vector F(nb_dof_rhs * N);
  for (size_type i = 0; i < nb_dof_rhs; ++i)
      gmm::copy(sol_f(mf_rhs.point_of_dof(i)),
                gmm::sub_vector(F, gmm::sub_interval(i*N, N)));
  
  // Volumic source term brick.
  getfem::mdbrick_source_term<> VOL_F(INCOMP, mf_rhs, F);
  
  // !!!NEW defining the left Dirichlet condition value.
  gmm::clear(F);
  getfem::mdbrick_Dirichlet<> left_bndy(VOL_F, LEFT_BOUNDARY_NUM);
  left_bndy.rhs().set(mf_rhs,F);

  // Defining the Dirichlet condition value.
  gmm::clear(F);
  getfem::mdbrick_Dirichlet<> final_model(left_bndy, DIRICHLET_BOUNDARY_NUM);
  final_model.rhs().set(mf_rhs,F);

  // Generic solve.
  cout << "Number of variables : " << final_model.nb_dof() << endl;
  getfem::standard_model_state MS(final_model);
  gmm::iteration iter(residual, 1, 40000);
  getfem::standard_solve(MS, final_model, iter);

  // Solution extraction
  gmm::copy(ELAS.get_solution(MS), U);
  
  return (iter.converged());
}

/**************************************************************************/
/*  main program.                                                         */
/**************************************************************************/

int main(int argc, char *argv[]) {

  DAL_SET_EXCEPTION_DEBUG; // Exceptions make a memory fault, to debug.
  FE_ENABLE_EXCEPT;        // Enable floating point exception for Nan.

  try {    
    stokes_problem p;
    p.PARAM.read_command_line(argc, argv);
    p.init();
    // p.mesh.write_to_file(p.datafilename + ".mesh");
    plain_vector U(p.mf_u.nb_dof());
    if (!p.solve(U)) DAL_THROW(dal::failure_error,"Solve has failed");
    if (p.PARAM.int_value("VTK_EXPORT")) {
      cout << "export to " << p.datafilename + ".vtk" << "..\n";
      getfem::vtk_export exp(p.datafilename + ".vtk",
                             p.PARAM.int_value("VTK_EXPORT")==1);
      exp.exporting(p.mf_u); 
      exp.write_point_data(p.mf_u, U, "stokes_displacement");
      cout << "export done, you can view the data file with (for example)\n"
        "mayavi -d " << p.datafilename << ".vtk -f ExtractVectorNorm -f "
        "WarpVector -m BandedSurfaceMap -m Outline\n";

      std::string outName = p.datafilename + ".out";
      std::ofstream os(outName.c_str());
      size_type dim = p.mf_u.get_qdim();
      cout << "Writing solution " << dim << ", for " << p.mf_u.nb_dof() << " 
dof\n";
      for (size_type i = 0; i < p.mf_u.nb_dof(); i += dim) {
        const base_node &node = p.mf_u.point_of_dof(i);
        for(size_type j = 0; j < node.size(); j++) 
          os << node[j] << " ";
        for(size_type j = 0; j < dim; j++) 
          os << U[i + j] << " ";
        os << "\n";
      }
      
    }
  }
  DAL_STANDARD_CATCH_ERROR;

  return 0; 
}