Hi Daniel,

are you by chance running 7.1.0 or an older version of deal.II? The
example compiles correctly with a recent development version for me.

On Wed, May 30, 2012 at 5:38 PM, Daniel Brauss <[email protected]> wrote:
> Dear all,
>
> I am trying to use extract_boundary_dofs for a  parallel distributed
> triangulation
> similar to step 40.  I will be setting up a constraint for the pressure term
> of the vector
> finite element being used ( similar to step 11 ) to obtain a unique pressure
> solution.
> I will be using extract_boundary_dofs to obtain the pressure dofs on the
> boundary
> of the domain where x=0 for a particular mpi process.  I have tried to
> simplify the
> code down to this one step on process 0 and have attached it below.
> I am receiving the following error during make
>
> aubkdb@dmc:ZeroBoundaryConstraint> make
> ==============debug========= ZeroBoundaryConstraintVersion2.cc  ->
> ZeroBoundaryConstraintVersion2.g.o
> ZeroBoundaryConstraintVersion2.cc: In member function ‘void
> MHD::MHDProblem::setup_system()’:
> ZeroBoundaryConstraintVersion2.cc:261: error: no matching function for call
> to ‘extract_boundary_dofs(dealii::DoFHandler<3, 3>&, std::vector<bool,
> std::allocator<bool> >&, dealii::IndexSet&, std::set<unsigned char,
> std::less<unsigned char>, std::allocator<unsigned char> >&)’
> ZeroBoundaryConstraintVersion2.cc:266: error: no matching function for call
> to ‘extract_boundary_dofs(dealii::DoFHandler<3, 3>&, std::vector<bool,
> std::allocator<bool> >&, dealii::IndexSet&, std::set<unsigned char,
> std::less<unsigned char>, std::allocator<unsigned char> >&)’
> make: *** [ZeroBoundaryConstraintVersion2.g.o] Error 1
> aubkdb@dmc:ZeroBoundaryConstraint>
>
>
>
> Thanks,
> Dan
>
>
>
>
> /* Include Files -----------------------------------------------------*/
>
> #include <base/quadrature_lib.h>
> #include <base/function.h>
> #include <base/timer.h>
> #include <base/utilities.h>               /* location for mpi
> queries              */
> #include <base/index_set.h>               /* dof indices working on
> processor */
> #include <base/conditional_ostream.h>     /* output from processor zero
> only       */
>
> #include <grid/grid_generator.h>
> #include <grid/tria_accessor.h>
> #include <grid/tria_iterator.h>
> #include <grid/grid_out.h>
> #include <grid/grid_tools.h>
>
> #include <distributed/tria.h>             /* parallel distributed
> Triangulation    */
> #include <distributed/grid_refinement.h>  /* parallel distributed grid
> refinement  */
>
> #include <dofs/dof_handler.h>
> #include <dofs/dof_accessor.h>
> #include <dofs/dof_tools.h>
>
> #include <fe/fe_q.h>
> #include <fe/fe_values.h>
> #include <fe/fe_system.h>
> #include <fe/fe_nedelec.h>
>
> #include <lac/vector.h>
> #include <lac/full_matrix.h>
> #include <lac/solver_cg.h>
> #include <lac/constraint_matrix.h>
> #include <lac/compressed_simple_sparsity_pattern.h>
> #include <lac/sparsity_tools.h>            /* parallel distributed sparsity
> pattern */
> #include <lac/petsc_parallel_sparse_matrix.h>
> #include <lac/petsc_parallel_vector.h>
> #include <lac/petsc_solver.h>
> #include <lac/petsc_precondition.h>
>
> #include <numerics/vectors.h>
> #include <numerics/data_out.h>
> #include <numerics/error_estimator.h>
>
> #include <fstream>
> #include <iostream>
>
> namespace MHD
> {
>
>   using namespace dealii;
>
>   class MHDProblem
>   {
>     public:
>       MHDProblem ();
>       ~MHDProblem ();
>
>       void run ();
>
>     private:
>       void setup_system ();
>
>       MPI_Comm mpi_communicator;
>
>       parallel::distributed::Triangulation<3>  triangulation;
>
>       DoFHandler<3>                            dof_handler;
>       FESystem<3>                              fe;
>
>       IndexSet                                 locally_owned_dofs;
>       IndexSet                                 locally_relevant_dofs;
>
>       ConstraintMatrix                         constraints;
>       // constraints to set one pressure and one potential basis function to
> zero
>       // also constraints for the dirichlet boundary conditions
>
>       PETScWrappers::MPI::SparseMatrix         system_matrix;
>       PETScWrappers::MPI::Vector               locally_relevant_solution;
>       PETScWrappers::MPI::Vector               system_rhs;
>
>       ConditionalOStream                       pcout;
>
>       const unsigned int                       n_mpi_processes;
>       const unsigned int
> this_mpi_process;
>
>   };
>
>
>
>
>
> /* MHD PROBLEM
> ********************************************************************/
> /**********************************************************************************/
> /**********************************************************************************/
>
> /* MHD Problem Constructor
> -------------------------------------------------*/
> /*--------------------------------------------------------------------------*/
>
>
>   MHDProblem::MHDProblem ()
>               :
>               mpi_communicator (MPI_COMM_WORLD),
>               triangulation (mpi_communicator,
>                              Triangulation<3>::MeshSmoothing
>                              (Triangulation<3>::smoothing_on_refinement |
>                               Triangulation<3>::smoothing_on_coarsening)),
>               dof_handler (triangulation),
>               fe(FE_Q<3>(2), 3,         // velocity
>                  FE_Q<3>(1), 1,         // pressure
>                  FE_Nedelec<3>(1), 1,   // current density
>                  FE_Q<3>(2), 1),        // electric potential
>               pcout (std::cout,
>                      (Utilities::System::
>                       get_this_mpi_process(mpi_communicator)
>                       == 0)),
>               n_mpi_processes
> (Utilities::MPI::n_mpi_processes(mpi_communicator)),
>               this_mpi_process
> (Utilities::MPI::this_mpi_process(mpi_communicator))
>   {}
>
>
>
> /* MHD Problem Destructor
> -------------------------------------------------*/
> /*--------------------------------------------------------------------------*/
>
>
>   MHDProblem::~MHDProblem ()
>   {
>     dof_handler.clear();
>   }
>
>
>
> /* MHD Problem System Setup
> ------------------------------------------------*/
> /*--------------------------------------------------------------------------*/
>
>
>   void MHDProblem::setup_system ()
>   {
>
>     dof_handler.distribute_dofs (fe);
>
>     pcout << "Distributed Degrees of Freedom (dofs)"
>           << std::endl;
>
>     locally_owned_dofs = dof_handler.locally_owned_dofs();
>
>     DoFTools::extract_locally_relevant_dofs (dof_handler,
>                                              locally_relevant_dofs);
>
>     locally_relevant_solution.reinit (mpi_communicator,
>                                       locally_owned_dofs,
>                                       locally_relevant_dofs);
>
>     locally_relevant_solution = 0;
>     system_rhs.reinit (mpi_communicator,
>                        dof_handler.n_dofs(),
>                        dof_handler.n_locally_owned_dofs());
>
>     system_rhs = 0;
>
>
>     /*-- Setting flags on boundary for dirichlet part and zero constraint
> -*/
>
>     unsigned char bD = '1';            // Dirichlet boundary flag
>     unsigned char x_equal_zero = '2';       // x component equal zero
> boundary flag
>
>     for (Triangulation<3>::active_cell_iterator
>             cell = dof_handler.begin_active();
>             cell != dof_handler.end(); ++cell)
>       if (cell->subdomain_id() == triangulation.locally_owned_subdomain())
>       {
>         for (unsigned int f=0; f < GeometryInfo<3>::faces_per_cell; ++f)
>         {
>           if (cell->face(f)->boundary_indicator() == 0)  // checking that
> face is on boundary
>           {
>             if (cell->face(f)->center()[0] == 0.0)       // x-component is 0
>               cell->face(f)->set_all_boundary_indicators(x_equal_zero);
>             if (cell->face(f)->center()[0] == 1.0)       // x-component is
> 1.0
>               cell->face(f)->set_all_boundary_indicators(bD);
>             if (cell->face(f)->center()[1] == 0.0)       // y-component is 0
>               cell->face(f)->set_all_boundary_indicators(bD);
>             if (cell->face(f)->center()[1] == 1.0)       // y-component is
> 1.0
>               cell->face(f)->set_all_boundary_indicators(bD);
>             if (cell->face(f)->center()[2] == 0.0)       // z-component is 0
>               cell->face(f)->set_all_boundary_indicators(bD);
>             if (cell->face(f)->center()[2] == 1.0)       // z-component is
> 1.0
>               cell->face(f)->set_all_boundary_indicators(bD);
>           }
>         }
>       }
>
>
>     pcout << "Set Boundary Flags " << std::endl;
>
>
>     if (this_mpi_process == 0)
>     {
>
>       std::set<unsigned char> boundary_indicators_x_equal_zero;
>       boundary_indicators_x_equal_zero.insert(x_equal_zero);
>       // declaring the boundary indicator(s) we will pass in to extract the
>       // basis functions having the corresponding boundary indicator flag
> (x_equal_zero)
>       // on the locally owned cells.  We will also be doing this for the
> pressure
>       // and electric potential using the masks show below
>
>
>       std::vector<bool> press_select (3*2 + 2, false);
>       press_select[3] = true;
>       // mask to select pressure basis function from current-density,
> potential, velocity basis fcns
>
>
>       IndexSet          boundary_dofs_pressure;
>       // setting up the index set that will be returned from
> extract_boundary_dofs
>       // giving the dofs that have a boundary indicator x_equal_zero and are
>       // pressure basis functions
>       //---- Not sure if the size of this set has to be specified and what
> the size
>       //     should be if it is only the locally owned cells -------
>
>       DoFTools::extract_boundary_dofs (dof_handler, press_select,
>                                        boundary_dofs_pressure,
>                                        boundary_indicators_x_equal_zero);
>       // getting the pressure basis functions that are on the boundary at x
> = 0
>
>     }
>
>   }
>
>
>   void MHDProblem::run ()
>   {
>
>     const unsigned int n_cycles = 1;
>
>     pcout << "this is process " << this_mpi_process
>           << std::endl << "the number of processes are "
>           << n_mpi_processes
>           << std::endl;
>
>     for (unsigned int cycle=0; cycle < n_cycles; ++cycle)
>     {
>       pcout << "Cycle " << cycle << ';' << std::endl;
>
>       if (cycle == 0)
>       {
>         const Point<3> bottom_left = Point<3> (0.0,0.0,0.0);
>         // bottom left corner of domain
>         const Point<3> top_right = Point<3> (1.0, 1.0, 1.0);
>         // top right corner of domain
>         std::vector<unsigned int> subdivisions (3,1);
>
>         subdivisions[0] = 2;  // divide domain in x-direction into two
>         subdivisions[1] = 2;  // divide domain in y-direction into two
>         subdivisions[2] = 2;  // divide domain in in z-direction
>         GridGenerator::subdivided_hyper_rectangle (triangulation,
>                                                  subdivisions,
>                                                  bottom_left,
>                                                  top_right);
>
>         triangulation.refine_global (0);
>
>         pcout << "Completed Grid Generation" << std::endl;
>
>         //std::ofstream output_file2("mesh");
>         char output_file2[]="mesh";
>         triangulation.write_mesh_vtk(output_file2);
>       }
>
>       setup_system ();
>
>
>       pcout << "  Number of active cells:"
>             << triangulation.n_global_active_cells()
>             << std::endl
>             << "  Number of degrees of freedom: "
>             << dof_handler.n_dofs()
>             << std::endl;
>
>     }
>
>     if (this_mpi_process == 0)
>     {
>       std::ofstream output_file("grid_out_p0.msh");
>       GridOut grid_out;
>       grid_out.write_msh (triangulation, output_file);
>     }
>
>   }
>
>
> }
>
> int main (int argc, char *argv[])
> {
>   try
>   {
>     using namespace dealii;
>     using namespace MHD;
>
>     PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
>     deallog.depth_console (0);
>
>     {
>       MHDProblem mhd_problem_3d;
>       mhd_problem_3d.run ();
>     }
>
>     PetscFinalize();
>
>   }
>   catch (std::exception &exc)
>   {
>     std::cerr << std::endl << std::endl
>               << "----------------------------------------"
>               << std::endl;
>     std::cerr << "Exception on processing: " << std::endl
>               << exc.what() << std::endl
>               << "Aborting!" << std::endl
>               << "----------------------------------------"
>               << std::endl;
>
>     return 1;
>   }
>   catch (...)
>   {
>     std::cerr << std::endl << std::endl
>               << "----------------------------------------"
>               << std::endl;
>     std::cerr << "Unknown exception!" << std::endl
>               << "Aborting!" << std::endl
>               << "----------------------------------------"
>               << std::endl;
>     return 1;
>   }
>
>   return 0;
> }
>
>
>
>
>
>
>
> _______________________________________________
> dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
>



-- 
Timo Heister
http://www.math.tamu.edu/~heister/
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii

Reply via email to