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
