#include <fstream>
#include <iostream>
-// #define DEBUG_OUTPUT_VTK
+//#define DEBUG_OUTPUT_VTK
+unsigned int counter=0;
const double eps = 1e-10;
void test2cells(const unsigned int p1=2,
const unsigned int p2=1)
{
+ Assert(dim>1,ExcInternalError());
Triangulation<dim> triangulation;
{
- Triangulation<dim> triangulationL;
- Triangulation<dim> triangulationR;
- GridGenerator::hyper_cube (triangulationL, -1,0); //create a square [-1,0]^d domain
- GridGenerator::hyper_cube (triangulationR, -1,0); //create a square [-1,0]^d domain
- Point<dim> shift_vector;
- shift_vector[0] = 1.0;
- GridTools::shift(shift_vector,triangulationR);
- GridGenerator::merge_triangulations (triangulationL, triangulationR, triangulation);
+ Point<dim> p1,p2;
+ for (unsigned int d=0; d<dim; d++)
+ p1[d] = -1;
+ p2[0] = 1.0;
+ std::vector<unsigned int> repetitoins(dim,1);
+ repetitoins[0]=2;
+ GridGenerator::subdivided_hyper_rectangle(triangulation,
+ repetitoins,
+ p1,
+ p2);
}
hp::DoFHandler<dim> dof_handler(triangulation);
#ifdef DEBUG_OUTPUT_VTK
// output to check if all is good:
+ counter++;
std::vector<Vector<double>> shape_functions;
std::vector<std::string> names;
for (unsigned int s=0; s < dof_handler.n_dofs(); s++)
for (unsigned int i = 0; i < shape_functions.size(); i++)
data_out.add_data_vector (shape_functions[i], names[i]);
- data_out.build_patches(15);
-
- std::string filename = "2cell-shape_functions_p1="+std::to_string(p1)
- +"_p2="+std::to_string(p2)
- +"_"+dealii::Utilities::int_to_string(dim)+"D.vtu";
+ data_out.build_patches(0);
+ std::string filename = "shape_functions_"+dealii::Utilities::int_to_string(counter,1)
+ +"_"+dealii::Utilities::int_to_string(dim)+"D.vtu";
std::ofstream output (filename.c_str ());
data_out.write_vtu (output);
#endif
std::vector<unsigned int> local_face_dof_indices;
std::vector<Vector<double> > values;
for (typename hp::DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active();
+ cell = dof_handler.begin_active();
cell != dof_handler.end(); ++cell)
- {
- const unsigned int fe_index = cell->active_fe_index();
- local_face_dof_indices.resize(fe_collection[fe_index].dofs_per_face);
- for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
- if (!cell->at_boundary(f)) // that's enough for our simple mesh
- {
- deallog << "cell="<<cell<<" face="<<f<<std::endl;
- fe_face_values_hp.reinit(cell,f);
- const FEFaceValues<dim> &fe_face_values = fe_face_values_hp.get_present_fe_values();
-
- const unsigned int n_q_points = fe_face_values.n_quadrature_points;
- values.resize (n_q_points,
- Vector<double>(2));
-
- fe_face_values.get_function_values (solution,values);
-
- const std::vector<dealii::Point<dim> > &q_points = fe_face_values.get_quadrature_points();
-
- for (unsigned int q = 0; q < n_q_points; q++)
- deallog << "u["<<q_points[q]<<"]={"<<values[q][0]<<","<<values[q][1]<<"}"<<std::endl;
- }
- }
+ {
+ const unsigned int fe_index = cell->active_fe_index();
+ local_face_dof_indices.resize(fe_collection[fe_index].dofs_per_face);
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (!cell->at_boundary(f)) // that's enough for our simple mesh
+ {
+ deallog << "cell="<<cell<<" face="<<f<<std::endl;
+ fe_face_values_hp.reinit(cell,f);
+ const FEFaceValues<dim> &fe_face_values = fe_face_values_hp.get_present_fe_values();
+
+ const unsigned int n_q_points = fe_face_values.n_quadrature_points;
+ values.resize (n_q_points,
+ Vector<double>(2));
+
+ fe_face_values.get_function_values (solution,values);
+
+ const std::vector<dealii::Point<dim> > &q_points = fe_face_values.get_quadrature_points();
+
+ for (unsigned int q = 0; q < n_q_points; q++)
+ deallog << "u["<<q_points[q]<<"]={"<<values[q][0]<<","<<values[q][1]<<"}"<<std::endl;
+ }
+ }
dof_handler.clear();
}
try
- {
- test2cells<2>(1,2); // Q(1) x FE_Nothing vs Q(2) x Q(1) => common Q(1) x FE_Nothing
- test2cells<2>(2,1); // Q(2) x FE_Nothing vs Q(1) x Q(1) => common Q(1) x FE_Nothing
- test2cells<3>(1,2); // Q(1) x FE_Nothing vs Q(2) x Q(1) => common Q(1) x FE_Nothing
- test2cells<3>(2,1); // Q(2) x FE_Nothing vs Q(1) x Q(1) => common Q(1) x FE_Nothing
- }
+ {
+ test2cells<2>(1,2); // Q(1) x FE_Nothing vs Q(2) x Q(1) => common Q(1) x FE_Nothing
+ test2cells<2>(2,1); // Q(2) x FE_Nothing vs Q(1) x Q(1) => common Q(1) x FE_Nothing
+ test2cells<3>(1,2); // Q(1) x FE_Nothing vs Q(2) x Q(1) => common Q(1) x FE_Nothing
+ test2cells<3>(2,1); // Q(2) x FE_Nothing vs Q(1) x Q(1) => common Q(1) x FE_Nothing
+ }
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;
- }
+ {
+ 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;
- };
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
}