Hi Timo, I am running 7.1.0 (double-checked). I have run step 40 successfully on the cluster along with testing a couple other extract functions.
I just looked at the DoFTools Namespace Reference on the web and see that only one extract_boundary_dofs function is listed for 7.1.0 and only a std::vector<bool> argument is available (rather than IndexSet) I guess I will have get a development version. Thanks. Dan On Thu, May 31, 2012 at 11:28 AM, Timo Heister <[email protected]>wrote: > 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
