Thanks Wolfgang,
I did as suggested and figured out the error. Although, I have also solved
it, I would like to know why this was a problem. For the same, I have
attached the new simple version. In that code, there is a piece of sub code
added which is actually the content of boundary_info() function which
existed earlier.
for (typename Triangulation<dim>::active_cell_iterator
cell = tria.begin();
cell != tria.end();
++cell ){
for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell;
++f){
Point<dim> point1 = cell->face(f)->center();
std::vector<double> x(dim);
for(unsigned int i = 0; i < dim; ++i){
x[i] = point1(i);
}
if (fabs(point1(0) - 0.0) < 1e-12)
cell->face(f)->set_boundary_indicator(42);
else
cell->face(f)->set_boundary_indicator(0);
}
}
If I remove the else part of the if-condition, it works, else the
constraints become zero. This implies, I should not try to set boundary
indicator for anything which is not at the boundary, which is obvious. I
know it was a mistake from my side, but I would like to know why it was a
problem (what goes on in the background as per the implementation?)
Best
Deepak
On Tue, Sep 27, 2016 at 5:02 PM, Wolfgang Bangerth <[email protected]>
wrote:
> On 09/27/2016 08:52 AM, Deepak Gupta wrote:
>
>>
>> I created a simple example but it worked correctly and is attached (I
>> modified
>> a test code given earlier my Jean-Paul). The difference between this
>> simple
>> code and my code stated earlier is that between the two checks of hanging
>> nodes, the interpolate_boundary_values function is not called in
>> hp_check.cpp.
>> However, in the piece of code sent earlier, there is this line:
>>
>> VectorTools::interpolate_boundary_values(dof_handler,
>> 42,
>> boundary_v,
>> boundary_values);
>>
>> This is somehow causing the problem. If I comment this one, the hanging
>> constraints are correctly. I am trying to put this as well in my
>> simplified
>> example. In the meantime, may you be you have an idea why this is messing
>> up.
>> The online documentation is down, so could not check details of the
>> function.
>>
>
> Well, VectorTools::interpolate_boundary_values() doesn't touch the
> 'constraints' variable, so it really shouldn't have any influence.
>
> When trying to find issues such as this, I usually try to start with
> something that doesn't work, and remove stuff as far as I can. So, start
> with your non-working program and make it smaller and smaller until the
> issue goes away. The part that you removed last was then the offending
> statement. Put it back in, and that is your minimal testcase.
>
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> www: http://www.math.colostate.edu/~bangerth/
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/fo
> rum/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 [email protected].
> 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 [email protected].
For more options, visit https://groups.google.com/d/optout.
#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>
#include <math.h>
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] = 3;
Point<dim> pt_1, pt_2;
pt_2[0] = 3;
pt_2[1] = 1;
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;
constraints.clear();
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
for (typename Triangulation<dim>::active_cell_iterator
cell = tria.begin();
cell != tria.end();
++cell ){
for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f){
Point<dim> point1 = cell->face(f)->center();
std::vector<double> x(dim);
for(unsigned int i = 0; i < dim; ++i){
x[i] = point1(i);
}
if (fabs(point1(0) - 0.0) < 1e-12)
cell->face(f)->set_boundary_indicator(42);
else
cell->face(f)->set_boundary_indicator(0);
}
}
constraints.close();
// 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";
//Second round of checking hanging nodes
dof_handler.distribute_dofs(fe_collection);
// Build constraints
constraints.clear();
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
constraints.close();
// 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;
}
