Last active
September 23, 2016 06:14
-
-
Save cpraveen/6c50ca2594aa0ead9ea2e95dc1e4c6a9 to your computer and use it in GitHub Desktop.
Testing sphere surface quadrature
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#include <deal.II/base/quadrature_lib.h> | |
#include <deal.II/grid/tria.h> | |
#include <deal.II/grid/tria_iterator.h> | |
#include <deal.II/grid/tria_accessor.h> | |
#include <deal.II/grid/grid_generator.h> | |
#include <deal.II/grid/grid_in.h> | |
#include <deal.II/grid/grid_out.h> | |
#include <deal.II/grid/tria_boundary_lib.h> | |
#include <deal.II/grid/grid_generator.h> | |
#include <deal.II/grid/grid_tools.h> | |
#include <deal.II/grid/manifold_lib.h> | |
#include <deal.II/dofs/dof_handler.h> | |
#include <deal.II/dofs/dof_accessor.h> | |
#include <deal.II/dofs/dof_tools.h> | |
#include <deal.II/fe/fe_q.h> | |
#include <deal.II/fe/fe_dgq.h> | |
#include <deal.II/fe/fe_values.h> | |
#include <deal.II/fe/mapping_q.h> | |
#include <deal.II/fe/mapping_q_eulerian.h> | |
#include <deal.II/fe/mapping_fe_field.h> | |
#include <deal.II/fe/fe_system.h> | |
#include <deal.II/numerics/vector_tools.h> | |
#include <cmath> | |
#include <iostream> | |
#include <fstream> | |
#include <string> | |
using namespace dealii; | |
int main(int argc, char ** argv) | |
{ | |
const unsigned int quad_order = 35; | |
const unsigned int quad_max = 50; | |
const double exact_area=4.*numbers::PI; | |
const unsigned int ref_max=2; | |
const unsigned int degree_max=14; | |
const double tol=1e-10; | |
std::cout<<"Testing the surface computation of a sphere"<<std::endl; | |
std::cout<<"exact area : "<<exact_area<<std::endl; | |
for(unsigned int i=1; i<=degree_max; ++i) | |
{ | |
std::cout<<"deg=" << i<<" "; | |
for(unsigned int j=i+1; j<=i+1; ++j) | |
{ | |
double surface_area=0.; | |
Triangulation<2,3> tria; | |
Point<3> center; | |
GridGenerator::hyper_sphere ( tria,center,1. ); | |
SphericalManifold<2, 3> manifold(center); | |
// Manifold<2,3> | |
// manifold = std::shared_ptr<SphericalManifold<2, 3> > (new SphericalManifold<2, 3> (pippo)); | |
tria.set_all_manifold_ids(0); | |
tria.set_manifold(0, manifold); | |
tria.refine_global(3); | |
FE_DGQArbitraryNodes<2,3> fe(QGaussLobatto<1>(i+1)); | |
FESystem<2, 3> fee(FE_Q<2, 3> (i), 3); | |
DoFHandler<2,3> dh(tria); | |
DoFHandler<2,3> dhh(tria); | |
dh.distribute_dofs(fe); | |
dhh.distribute_dofs(fee); | |
std::cout<<"ndofs=" << dh.n_dofs()<<" "; | |
Vector<double> euler_vec(dhh.n_dofs()); | |
VectorTools::get_position_vector(dhh,euler_vec); | |
// MappingFEField<2,3> mapping(dhh,euler_vec); | |
MappingQ<2,3> mapping(i); | |
QGaussLobatto<2> quad(j); | |
QGaussLobatto<1> face_quad(j); | |
const unsigned int n_q_points = quad.size(); | |
FEValues<2,3> fe_values (mapping, fe, quad, | |
update_values | | |
update_quadrature_points | | |
update_JxW_values); | |
FEFaceValues<2,3> fe_face_values (mapping, fe, face_quad, | |
update_values | | |
update_quadrature_points | | |
update_normal_vectors | | |
update_JxW_values); | |
FEFaceValues<2,3> fe_face_values_nbr (mapping, fe, face_quad, | |
update_values | | |
update_quadrature_points | | |
update_normal_vectors | | |
update_JxW_values); | |
std::vector<Point<3> > support_points(dh.n_dofs()); | |
DoFTools::map_dofs_to_support_points<2, 3>(mapping, | |
dh, support_points); | |
for(unsigned int i=0;i<support_points.size();++i) | |
if(std::fabs(support_points[i].square()-1.)>tol) | |
{ | |
std::cout<<"ERROR!!! : "<<std::fabs(support_points[i].square()-1.)<<std::endl; | |
break; | |
} | |
// for(unsigned int i=0;i<support_points.size();++i) | |
// if(std::fabs(support_points[i].square()-1.)<tol) | |
// { | |
// std::cout<<"OK!!! : "<<std::fabs(support_points[i].square()-1.)<<std::endl; | |
// } | |
double max_dw= 0, max_dn = 0; | |
for (typename DoFHandler<2,3>::active_cell_iterator | |
cell = dh.begin_active(), | |
endc = dh.end(); | |
cell!=endc; ++cell) | |
{ | |
fe_values.reinit(cell); | |
for (unsigned int q_point=0; q_point<n_q_points; ++q_point) | |
{ | |
surface_area += fe_values.JxW(q_point); | |
} | |
for(unsigned int f=0; f<GeometryInfo<2>::faces_per_cell; ++f) | |
{ | |
fe_face_values.reinit(cell, f); | |
fe_face_values_nbr.reinit (cell->neighbor(f), cell->neighbor_of_neighbor(f)); | |
for (unsigned int q_point=0; q_point<face_quad.size(); ++q_point) | |
{ | |
double dw = std::fabs(fe_face_values.JxW(q_point) - fe_face_values_nbr.JxW(q_point)); | |
dw /= fe_face_values.JxW(q_point); | |
max_dw = std::max(max_dw, dw); | |
Tensor<1,3> dn = fe_face_values.normal_vector(q_point) + | |
fe_face_values_nbr.normal_vector(q_point); | |
max_dn = std::max(max_dn, dn.norm()); | |
} | |
} | |
} | |
// std::cout<<surface_area<<" "<<exact_area<<std::endl; | |
double error = std::fabs(surface_area-exact_area)/exact_area; | |
// std::cout<<"degree "<<i<<" error "<<error<<std::endl; | |
std::cout<<", area error =" << error<<std::endl; | |
std::cout << "max dweight =" << max_dw << " max dnormal=" << max_dn << std::endl; | |
std::string output_filename; | |
output_filename="sphere_order_"+Utilities::int_to_string(i)+"_refinement_"+Utilities::int_to_string(j)+".vtu"; | |
std::ofstream output_file(output_filename); | |
GridOut().write_vtu (tria, output_file); | |
tria.set_manifold(0); | |
} | |
std::cout<<std::endl; | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Since this issue seems to disappear when changing
w
into1-w
, as you indicate in one of the issues, my guess is that there is an internal problem in the computation of the "intermediate support points" of the finite elements.I'm looking into this...