Hi, I'm working on a flux-pressure model using a Raviart-Thomas/DGQ FESystem. I'm using Dirichlet conditions on the RT (flux) elements on part of the domain and on the DGQ (pressure) elements on another part of the domain using the boundary_id's. I'm using VectorTools::interpolate_boundary_values to add the pressure BCs, but it doesn't appear to be adding rows to the ConstraintMatrix.
The enclosed minimal code shows that interpolate_boundary_values is computing the BC value function, but it also shows that there is no change in the ConstraintMatrix after it is called. I can't figure out what logic in do_interpolate_boundary_values is keeping it from adding the appropriate entries to the ConstraintMatrix. I would appreciate if anyone can find my error here or explain why I'm not getting the constraints I'm expecting. Thanks, Jonathan -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
//Jonathan L. Matthews San Diego State University, 2018
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_raviart_thomas.h>
#include <deal.II/numerics/vector_tools.h>
#include <iostream>
namespace LA
{
using namespace dealii::LinearAlgebraPETSc;
}
using std::cout;
namespace ResFlow
{
using namespace dealii;
template<int dim>
class PressureBoundaryValues : public Function<dim>
{
public:
PressureBoundaryValues () :
Function<dim> (dim+1)
{
}
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
virtual
~PressureBoundaryValues ()
{
}
;
};
template <int dim>
void PressureBoundaryValues<dim>::vector_value (const Point<dim> &p, Vector<double> &values) const
{
for (unsigned int c=0; c < this->n_components; ++c)
values(c) = c == dim ? 150. : 0.;
std::cout << "PressureBoundaryValues used to compute vector_value (for interpolate_boundary_values.)" << std::endl;
values.print(std::cout);
}
template<int dim>
class FluxBoundaryValues : public Function<dim>
{
public:
FluxBoundaryValues (Vector<double> flux) :
Function<dim> (1),
flux_density(flux)
{
}
FluxBoundaryValues () :
Function<dim> (1)
{
}
virtual void
vector_value (const Point<dim> &p, Vector<double> &flux) const;
void set_flux_density(Vector<double> flux)
{
flux_density = flux;
}
virtual
~FluxBoundaryValues ()
{
}
private:
Vector<double> flux_density;
};
template<int dim>
void
FluxBoundaryValues<dim>::vector_value (const Point<dim> & p, Vector<double> & flux) const
{
flux = flux_density;
}
template<int dim>
class ResFlowProblem
{
public:
ResFlowProblem (
const unsigned int degree,
Triangulation<dim> & triangulation);
~ResFlowProblem ();
void make_grid (Point<dim> p1, Point<dim> p2, unsigned int nx,
unsigned int ny, unsigned int nz);
void run ();
void setup_system ();
void set_inlet_flux_id(types::boundary_id id_no)
{
inlet_flux_boundary_id = id_no;
}
void set_pressure_bc_id(types::boundary_id id_no)
{
pressure_boundary_id = id_no;
}
void set_inlet_flux(Vector<double> flux);
const unsigned int degree;
FESystem<dim> fe;
private:
void make_grid ();
SmartPointer<Triangulation<dim> > triangulation;
DoFHandler<dim> dof_handler;
std::vector<IndexSet> owned_partitioning;
std::vector<IndexSet> relevant_partitioning;
ConstraintMatrix constraints;
types::boundary_id inlet_flux_boundary_id;
types::boundary_id pressure_boundary_id;
FluxBoundaryValues<dim> influx;
PressureBoundaryValues<dim> pressure_bc;
};
template<int dim>
ResFlowProblem<dim>::ResFlowProblem (
const unsigned int degree,
Triangulation<dim> & grid) :
degree (degree),
fe (FE_RaviartThomas<dim> (degree), 1, FE_DGQ<dim> (degree), 1),
triangulation (&grid),
dof_handler (*triangulation),
inlet_flux_boundary_id (0),
pressure_boundary_id(0)
{
}
template<int dim>
ResFlowProblem<dim>::~ResFlowProblem ()
{
dof_handler.clear ();
}
template<int dim>
void
ResFlowProblem<dim>::make_grid ()
{
// default make_grid function for testing
unsigned int nx = 8, ny = 8, nz = 8;
bool colorize = true;
std::vector<unsigned int> elem_dim;
if (dim == 2)
elem_dim =
{ nx,nz};
else
elem_dim =
{ nx,ny,nz};
Point<dim> p1, p2;
if (dim == 2)
{
p1 = Point<dim> (0., -153750.);
p2 = Point<dim> (50000., -151750.);
}
else
{
p1 = Point<dim> (0., 0., -153750.);
p2 = Point<dim> (50000., 50000., -151750.);
}
GridGenerator::subdivided_hyper_rectangle (*triangulation, elem_dim, p1,
p2, colorize);
} // end make_grid()
template<int dim>
void
ResFlowProblem<dim>::setup_system ()
{
dof_handler.distribute_dofs (fe);
std::vector<unsigned int> resflow_sub_blocks (2, 0);
resflow_sub_blocks[1] = 1;
// Renumber to yield block structure
DoFRenumbering::block_wise (dof_handler);
std::vector<types::global_dof_index> dofs_per_block (2);
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block,
resflow_sub_blocks);
const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1];
cout << " Number of degrees of freedom: " << dof_handler.n_dofs ()
<< " (" << n_u << '+' << n_p << ')' << std::endl;
IndexSet locally_relevant_dofs;
DoFTools::extract_locally_relevant_dofs (dof_handler,
locally_relevant_dofs);
{
std::cout << "\nConstraints before clearing: " << std::endl;
constraints.print(std::cout);
constraints.clear();
std::cout << "\nConstraints after clearing: " << std::endl;
constraints.print(std::cout);
constraints.reinit (locally_relevant_dofs);
std::cout << "\nConstraints after reinit: " << std::endl;
constraints.print(std::cout);
Assert(inlet_flux_boundary_id!=0,
ExcMessage("inlet_flux_boundary_id must be initialized and may not be set to zero. \n\nSet the boundary ID before calling ResFlowProblem<dim>::setup_system."));
VectorTools::project_boundary_values_div_conforming (dof_handler, 0,
influx, inlet_flux_boundary_id,
constraints);
std::cout << "\nConstraints after project_boundary_values_div_conforming: " << std::endl;
constraints.print(std::cout);
Assert(pressure_boundary_id!=0,
ExcMessage("pressure_boundary_id must be initialized and may not be set to zero. \n\nSet the boundary ID before calling ResFlowProblem<dim>::setup_system."));
FEValuesExtractors::Scalar pressure(dim);
ComponentMask pressure_mask = fe.component_mask (pressure);
VectorTools::interpolate_boundary_values(dof_handler,
pressure_boundary_id,
pressure_bc,
constraints,
pressure_mask);
std::cout << "\nConstraints after interpolate_boundary_values: " << std::endl;
constraints.print(std::cout);
constraints.close ();
std::cout << "\nConstraints after close: " << std::endl;
constraints.print(std::cout);
}
} // end setup_system
template<int dim>
void ResFlowProblem<dim>::set_inlet_flux(Vector<double> flux)
{
influx.set_flux_density(flux);
}
template<int dim>
void
ResFlowProblem<dim>::run ()
{
cout << "Running in serial." << std::endl;
make_grid ();
cout << "Grid made.\n" << std::endl;
set_inlet_flux_id(1);
Vector<double> inlet_flux(3);
inlet_flux(0) = 1.0;
set_inlet_flux(inlet_flux);
set_pressure_bc_id(3);
setup_system ();
cout << "\nSystem set up." << std::endl;
}
} // end namespace ResFlow
int
main (int argc, char *argv[])
{
const unsigned int dim = 2;
try
{
using namespace dealii;
using namespace ResFlow;
Assert(dim==2 || dim==3, ExcNotImplemented());
Triangulation<dim> grid;
ResFlowProblem<dim> resflow (0, grid);
resflow.run ();
}
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;
}
output
Description: Binary data
