From: bangerth Date: Mon, 20 Dec 2010 03:51:23 +0000 (+0000) Subject: Fix computation of face normal vectors in the codim=1 case. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6b10e6ff576189ae75bcb2cd4d56af620ff1343a;p=dealii-svn.git Fix computation of face normal vectors in the codim=1 case. git-svn-id: https://svn.dealii.org/trunk@23015 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/source/fe/mapping_q1.cc b/deal.II/source/fe/mapping_q1.cc index 200ed752ac..34f474b93b 100644 --- a/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/source/fe/mapping_q1.cc @@ -1016,7 +1016,8 @@ namespace internal Assert (JxW_values.size() == n_q_points, ExcDimensionMismatch(JxW_values.size(), n_q_points)); - // map the unit tangentials to the real cell + // map the unit tangentials to the + // real cell for (unsigned int d=0; d::faces_per_cell*d < @@ -1032,35 +1033,81 @@ namespace internal mapping_contravariant); } - switch (dim) - { - case 2: - { - for (unsigned int i=0; i cell_normal; + + data.contravariant[point] = transpose(data.contravariant[point]); + + // compute the normal + // vector to this cell + // and put it into the + // last row of + // data.contravariant + if ( (dim==1) && (spacedim==2) ) + cross_product(cell_normal, + -data.contravariant[point][0]); + else if ( (dim==2) && (spacedim==3) ) + cross_product(cell_normal, + data.contravariant[point][0], + data.contravariant[point][1]); + else + Assert (false, ExcNotImplemented()); + + // then compute the face + // normal from the face + // tangent and the cell + // normal: + if ( (dim==1) && (spacedim==2) ) + // need to think how to + // figure out whether we + // need to point to the + // left or right + Assert (false, ExcNotImplemented()) + else if ( (dim==2) && (spacedim==3) ) + cross_product (boundary_forms[point], + data.aux[0][point], cell_normal); + else + Assert (false, ExcNotImplemented()); + } + } + if (update_flags & (update_normal_vectors | update_JxW_values)) for (unsigned int i=0;i::subface_ratio( @@ -1070,7 +1117,7 @@ namespace internal } if (update_flags & update_normal_vectors) - normal_vectors[i] = boundary_forms[i] / f; + normal_vectors[i] = boundary_forms[i] / boundary_forms[i].norm(); } } } diff --git a/tests/codim_one/normal_vectors_01.cc b/tests/codim_one/normal_vectors_01.cc index 7b6a0160cf..84cda10e76 100644 --- a/tests/codim_one/normal_vectors_01.cc +++ b/tests/codim_one/normal_vectors_01.cc @@ -11,9 +11,8 @@ // //---------------------------- normal_vectors_01.cc --------------------------- -/* - Asking for normal vectors in the codim-1 case led to aborts. -*/ +// Asking for face normal vectors in the codim-1 case led to aborts. + #include "../tests.h" @@ -34,28 +33,11 @@ template void test () { - Triangulation tria; - std::map::cell_iterator, - typename Triangulation::face_iterator> surface_to_volume_mapping; - - HyperBallBoundary boundary_description; - Triangulation volume_mesh; - GridGenerator::half_hyper_ball(volume_mesh); - - volume_mesh.set_boundary (1, boundary_description); - volume_mesh.set_boundary (0, boundary_description); - volume_mesh.refine_global (1); - - static HyperBallBoundary surface_description; - tria.set_boundary (1, surface_description); - tria.set_boundary (0, surface_description); - - std::set boundary_ids; - boundary_ids.insert(0); + Triangulation volume_mesh; + GridGenerator::hyper_cube(volume_mesh); - surface_to_volume_mapping - = GridTools::extract_boundary_mesh (volume_mesh, tria, - boundary_ids); + Triangulation tria; + GridTools::extract_boundary_mesh (volume_mesh, tria); FE_Q fe (1); DoFHandler dh(tria); @@ -63,14 +45,27 @@ void test () dh.distribute_dofs (fe); - FEFaceValues fe_face_values (mapping, fe, QGauss(2), - update_values | update_quadrature_points | - update_normal_vectors | update_JxW_values); - - typename DoFHandler::active_cell_iterator - cell = dh.begin_active(); - - fe_face_values.reinit (cell,0); + FEFaceValues fe_face_values (mapping, fe, QTrapez(), + update_normal_vectors); + + for (typename DoFHandler::active_cell_iterator + cell = dh.begin_active(); + cell != dh.end(); ++cell) + { + deallog << "Face centered at " << cell->center() + << std::endl; + + for (unsigned int f=0; f::faces_per_cell; ++f) + { + deallog << " Edge centered at " << cell->face(f)->center() + << std::endl; + + fe_face_values.reinit (cell,f); + for (unsigned int q=0; q