Hi Deepak,
Ok, so I've created a minimal working example to demonstrate what you're
seeing here. I'm not sufficiently familiar with the p-refinement part of
the hp functionality/implementation to tell whether this is truly the
expected behaviour or not. I see that the examples and description
<https://www.dealii.org/8.4.1/doxygen/deal.II/group__constraints.html#ga3eaa31a679484e80c193e74e8a967dc8>
in the documentation mix h- and p-refinement, so are *perhaps* not
sufficiently clear to describe the expected outcome of this configuration.
Or maybe its a bug arising from the adjacent cells being on the same
refinement level (i.e. there are no actual hanging *nodes*). Perhaps
someone else could clarify this?
J-P
On Tuesday, August 23, 2016 at 9:45:41 AM UTC+2, Deepak Gupta wrote:
>
> Dear Daniel and others,
>
> Below is the part of my code that might explain better:
>
> *//Below is the constructor for FEM class (a part of my own library). The
> part of the code below in BOLD states that I assign different order shape
> functions alternatively.*
> template<int dim>
> FEM<dim>::FEM(
> Triangulation<dim> &obj_triangulation,
> hp::DoFHandler<dim> &dof_handler,
> DefineMesh<dim> &obj_mesh){
> this->mesh = &obj_mesh;
> this->dof_handler = &dof_handler; // for the state field on the
> analysis
> this->triangulation = &obj_triangulation; //for the state field on
> the analysis
>
> //For Lagrange element
> if (mesh->elementType == "FE_Q"){
> for (unsigned int degree = 1; degree <= mesh->max_el_order;
> ++degree){
> fe_collection.push_back(FESystem<dim>(FE_Q<dim>(degree), dim));
> }
>
> }
> //Quadrature collection for FE
> for (unsigned int qrule = 1; qrule <= mesh->max_el_order + 11;
> ++qrule){
> quadrature_collection.push_back(QGauss<dim>(qrule));
> }
> //initializing the type of element for each cell of analysis mesh
>
>
>
>
>
> * for (typename hp::DoFHandler<dim>::active_cell_iterator cell =
> dof_handler.begin_active(); cell != dof_handler.end(); ++cell){
> cell->set_active_fe_index(p_index); ++cell;
> cell->set_active_fe_index(p_index+2); }*
> }
>
>
> *//The setup() function is as follows, The line that prints the no. of
> constraints is put in BOLD in the code *template <int dim>
> void FEM<dim>::setup_system(){
> //FE mesh
> dof_handler->distribute_dofs(fe_collection);
> solution.reinit(dof_handler->n_dofs());
> system_rhs.reinit(dof_handler->n_dofs());
>
> analysis_density_handler->distribute_dofs(fe_analysis_density_collection);
> //Used to add density on every node
> hanging_node_constraints.clear();
> DoFTools::make_hanging_node_constraints(*dof_handler,
> hanging_node_constraints);
> hanging_node_constraints.close();
>
> *std::cout<<"No. of hanging node constraints :
> "<<hanging_node_constraints.n_constraints()<<std::endl;*
> DynamicSparsityPattern dsp (dof_handler->n_dofs());
> DoFTools::make_sparsity_pattern (*dof_handler,
> dsp,
> hanging_node_constraints,
> true);
> sparsity_pattern.copy_from(dsp);
> system_matrix.reinit(sparsity_pattern);
> }
>
>
> *OUTPUT: *No. of hanging node constraints : 0
>
>
> * ---------------------------------------------------------------------
> *Could
> someone let me know why the no. of constraints is zero. Since the p-order
> of neighboring elements differ, I would expect constraints here.
>
> Best regards
> Deepak
>
> On Fri, Aug 19, 2016 at 10:51 AM, Daniel Arndt <
> d.ar...@math.uni-goettingen.de> wrote:
>
>> Deepak,
>>
>> I am aware of the fact that ConstraintMatrix.n_constraints() gives the
>>> number of hanging nodes for h-refinement, but if I only use p-refinement,
>>> can it give the number of hanging support points occurring due to the
>>> difference in the order of bases between two adjacent elements? I am
>>> expecting this based on the description of the hp-refinement document.
>>>
>>> I have been trying this but always get 0 hanging support points,
>>> although I know that the some elements of the mesh use Q1, others Q2.
>>>
>> That sounds definitely strange. step-27 uses a hp::DoFHandler and
>> make_hanging_node_constraints. Do you see that constraints are created
>> there?
>> If you can't find a solution looking at that example program, can you
>> reduce your code to a minimal example showing that no constraints are
>> created for p-refinement?
>>
>> Best,
>> Daniel
>> --
>> 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.
>>
>
>
--
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.
##
# CMake script for the step-1 tutorial program:
##
# Set the name of the project and target:
SET(TARGET hp-pure_p_refinement)
# Declare all source files the target consists of:
SET(TARGET_SRC
${TARGET}.cc
# You can specify additional files here!
)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
FIND_PACKAGE(deal.II 8.4 QUIET
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/numerics/data_out.h>
#include <iostream>
#include <fstream>
using namespace dealii;
int main()
{
const int dim = 2;
Triangulation<dim> tria;
const hp::FECollection<dim> fe_collection (FE_Q<dim>(1), FE_Q<dim>(2));
hp::DoFHandler<dim> dof_handler(tria);
// Two cell domain
std::vector<unsigned int> divisions (dim,1);
divisions[0] = 2;
Point<dim> pt_1, pt_2;
for (unsigned int d=0; d<dim; ++d)
pt_2[d] = 1.0;
GridGenerator::subdivided_hyper_rectangle(tria, divisions, pt_1, pt_2);
// Set active FE index (p-refinement level; default index is 0)
dof_handler.begin_active()->set_active_fe_index(1);
dof_handler.distribute_dofs(fe_collection);
// Build constraints
ConstraintMatrix constraints;
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
constraints.clear();
// Print to screen
std::cout
<< "Number of cells: " << tria.n_active_cells() << "\n"
<< "Number of degrees of freedom: " << dof_handler.n_dofs() << "\n"
<< "Number of constraints: " << constraints.n_constraints() << "\n";
// Output p-refinement level
Vector<double> p_refinement (tria.n_active_cells());
unsigned int cell_no = 0;
for (typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active();
cell != dof_handler.end(); ++cell, ++cell_no)
{
p_refinement[cell_no] = cell->active_fe_index();
}
DataOut<dim,hp::DoFHandler<dim> > data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (p_refinement, "p_refinement_level");
data_out.build_patches ();
std::ostringstream filename;
filename << "solution.vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
return 0;
}
