From b220b4b078022abe3c699d6eed443416425d4634 Mon Sep 17 00:00:00 2001 From: guido Date: Tue, 6 Mar 2001 12:20:49 +0000 Subject: [PATCH] mapping merge git-svn-id: https://svn.dealii.org/trunk@4119 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/fe/Makefile.in | 7 + tests/fe/fe_data_test.cc | 78 ++++++ tests/fe/mapping.cc | 502 +++++++++++++++++++++++++++++++++++++++ tests/fe/shapes.cc | 241 +++++++++++++++++++ 4 files changed, 828 insertions(+) create mode 100644 tests/fe/fe_data_test.cc create mode 100644 tests/fe/mapping.cc create mode 100644 tests/fe/shapes.cc diff --git a/tests/fe/Makefile.in b/tests/fe/Makefile.in index b44d0abb69..7f60aa37bf 100644 --- a/tests/fe/Makefile.in +++ b/tests/fe/Makefile.in @@ -89,6 +89,7 @@ fe_data_test-o-files = $(fe_data_test-cc-files:.cc=.o) endif fe_data_test: $(fe_data_test-o-files) $(libraries) + @echo =====linking======= $< @$(CXX) $(flags) -o $@ $^ ############################################################ @@ -102,6 +103,7 @@ shapes-o-files = $(shapes-cc-files:.cc=.o) endif shapes: $(shapes-o-files) $(libraries) + @echo =====linking======= $< @$(CXX) $(flags) -o $@ $^ ############################################################ @@ -115,6 +117,7 @@ derivatives-o-files = $(derivatives-cc-files:.cc=.o) endif derivatives: $(derivatives-o-files) $(libraries) + @echo =====linking======= $< @$(CXX) $(flags) -o $@ $^ ############################################################ @@ -128,6 +131,7 @@ embedding-o-files = $(embedding-cc-files:.cc=.o) endif embedding: $(embedding-o-files) $(libraries) + @echo =====linking======= $< @$(CXX) $(flags) -o $@ $^ ############################################################ @@ -154,6 +158,7 @@ mapping-o-files = $(mapping-cc-files:.cc=.o) endif mapping: $(mapping-o-files) $(libraries) + @echo =====linking======= $< @$(CXX) $(flags) -o $@ $^ ############################################################ @@ -167,6 +172,7 @@ internals-o-files = $(internals-cc-files:.cc=.o) endif internals: $(internals-o-files) $(libraries) + @echo =====linking======= $< @$(CXX) $(flags) -o $@ $^ ############################################################ @@ -180,6 +186,7 @@ performance-o-files = $(performance-cc-files:.cc=.o) endif performance: $(performance-o-files) $(libraries) + @echo =====linking======= $< @$(CXX) $(flags) -o $@ $^ ############################################################ diff --git a/tests/fe/fe_data_test.cc b/tests/fe/fe_data_test.cc new file mode 100644 index 0000000000..cd7c3810d6 --- /dev/null +++ b/tests/fe/fe_data_test.cc @@ -0,0 +1,78 @@ +//---------------------------- fe_data_test.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- fe_data_test.cc --------------------------- + + +#include +#include + +#include +#include +#include +#include +#include + + +template +void test_fe_datas() +{ + vector *> fe_datas; + fe_datas.push_back(new FE_Q (1)); + fe_datas.push_back(new FE_Q (4)); + fe_datas.push_back(new FESystem(FE_Q (2), 2)); + fe_datas.push_back(new FESystem(FE_Q (1), 2, + FE_Q (2), 1)); + // for this an assertion thrown in + // FESystem::build_interface_constraints + // for the following calls for dim=3 +// fe_datas.push_back(new FESystem(FE_Q (3), 2)); +// fe_datas.push_back(new FESystem(FE_Q (1), 2, +// FE_Q (3), 1)); +// fe_datas.push_back(new FESystem(FE_Q (4), 2)); + + deallog << endl << "dim=" << dim << endl; + for (unsigned int n=0; n *fe_data=fe_datas[n]; + deallog << "fe_data[" << n <<"]:" << endl; + deallog << "dofs_per_vertex=" << fe_data->dofs_per_vertex << endl; + deallog << "dofs_per_line=" << fe_data->dofs_per_line << endl; + deallog << "dofs_per_quad=" << fe_data->dofs_per_quad << endl; + deallog << "dofs_per_hex=" << fe_data->dofs_per_hex << endl; + deallog << "first_line_index=" << fe_data->first_line_index << endl; + deallog << "first_quad_index=" << fe_data->first_quad_index << endl; + deallog << "first_hex_index=" << fe_data->first_hex_index << endl; + deallog << "first_face_line_index=" << fe_data->first_face_line_index << endl; + deallog << "first_face_quad_index=" << fe_data->first_face_quad_index << endl; + deallog << "dofs_per_face=" << fe_data->dofs_per_face << endl; + deallog << "dofs_per_cell=" << fe_data->dofs_per_cell << endl; + deallog << "components=" << fe_data->components << endl; + } + + // delete all FiniteElementDatas + for (unsigned int i=0; i(); + test_fe_datas<2>(); + test_fe_datas<3>(); +} + + + diff --git a/tests/fe/mapping.cc b/tests/fe/mapping.cc new file mode 100644 index 0000000000..ff693b26f6 --- /dev/null +++ b/tests/fe/mapping.cc @@ -0,0 +1,502 @@ +// $Id$ +// Copyright (C) 2001 Ralf Hartmann +// +// Show the shape functions implemented and computes the area of cells. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + +char fname[50]; + +template +inline void +plot_transformation(Mapping &mapping, + FiniteElement &fe, + DoFHandler::cell_iterator &cell, + const char* name) +{ + const unsigned int div = 7; + + QTrapez<1> q_trapez; + QIterated q(q_trapez, div); + FEValues fe_values(mapping, fe, q, + UpdateFlags(update_q_points + | update_JxW_values)); + + fe_values.reinit(cell); + const vector &JxW=fe_values.get_JxW_values(); + + ofstream gnuplot(name); + + unsigned int k=0; + for (unsigned int nz=0; nz<=((dim>2) ? div : 0); ++nz) + { + for (unsigned int ny=0; ny<=((dim>1) ? div : 0); ++ny) + { + for (unsigned int nx=0; nx<=div; ++nx) + { + gnuplot << fe_values.quadrature_point(k); + double J = JxW[k] / q.weight(k); + gnuplot << ' ' << J << endl; + ++k; + } + gnuplot << endl; + } + gnuplot << endl; + } +} + + + +template +inline void +plot_faces(Mapping &mapping, + FiniteElement &fe, + DoFHandler::cell_iterator &cell, + const char* name) +{ + ofstream gnuplot(name); + + QGauss4 q; + const unsigned int nq = (unsigned int) (.01 + pow(q.n_quadrature_points, 1./(dim-1))); + + FEFaceValues fe_values(mapping, fe, q, + UpdateFlags(update_q_points + | update_normal_vectors)); + + for (unsigned int face_nr=0; + face_nr < GeometryInfo::faces_per_cell; + ++ face_nr) + { + fe_values.reinit(cell, face_nr); + + const vector > &normals + =fe_values.get_normal_vectors(); + + unsigned int k=0; + for (unsigned int ny=0; ny<((dim>2) ? nq : 1); ++ny) + { + for (unsigned int nx=0; nx x = fe_values.quadrature_point(k); + Tensor<1,dim> n = normals[k]; + gnuplot << x << '\t' << n << endl; + ++k; + } + gnuplot << endl; + } + gnuplot << endl; + } +} + + + +template +inline void +plot_subfaces(Mapping &mapping, + FiniteElement &fe, + DoFHandler::cell_iterator &cell, + const char* name) +{ + ofstream gnuplot(name); + + QGauss4 q; + const unsigned int nq = (unsigned int) (.01 + pow(q.n_quadrature_points, 1./(dim-1))); + + FESubfaceValues fe_values(mapping, fe, q, + UpdateFlags(update_q_points + | update_normal_vectors)); + for (unsigned int face_nr=0; + face_nr < GeometryInfo::faces_per_cell; + ++ face_nr) + for (unsigned int sub_nr=0; + sub_nr < GeometryInfo::subfaces_per_face; + ++ sub_nr) + { + fe_values.reinit(cell, face_nr, sub_nr); + + const vector > &normals + =fe_values.get_normal_vectors(); + + unsigned int k=0; + for (unsigned int ny=0; ny<((dim>2) ? nq : 1); ++ny) + { + for (unsigned int nx=0; nx x = fe_values.quadrature_point(k); + Tensor<1,dim> n = normals[k]; + gnuplot << x << '\t' << n << endl; + ++k; + } + gnuplot << endl; + } + gnuplot << endl; + } +} + + + +template<> +inline void +plot_faces(Mapping<1>&, + FiniteElement<1>&, + DoFHandler<1>::cell_iterator&, + const char*) +{}; + + + +template<> +inline void +plot_subfaces(Mapping<1>&, + FiniteElement<1>&, + DoFHandler<1>::cell_iterator&, + const char*) +{}; + + + +template +inline void +compute_area(Mapping &mapping, + FiniteElement &fe, + DoFHandler::cell_iterator &cell) +{ + QGauss4 gauss4; + FEValues fe_values(mapping, fe, gauss4, + UpdateFlags(update_JxW_values)); + fe_values.reinit(cell); + const vector &JxW=fe_values.get_JxW_values(); + + double area=0; + for (unsigned int i=0; i *> tria_ptr; + vector *> boundary_ptr; + vector exact_areas; + + create_triangulations(tria_ptr, boundary_ptr, exact_areas); + Assert(show.size()==tria_ptr.size(), ExcInternalError()); + + FE_Q fe_q4(4); + + for (unsigned int i=0; i dof(*tria_ptr[i]); + dof.distribute_dofs(fe_q4); + DoFHandler::cell_iterator cell = dof.begin_active(); + + deallog << "Triangulation" << i << ":" << endl; + + deallog << "exact_area=" << exact_areas[i] << endl; + for (unsigned int j=0; j1) + { + ostrstream ost(st2, 99); + ost << "MappingFace" << dim << "d-" << i << '-' + << mapping_strings[j] << ".dat" << ends; + deallog << st2 << endl; + plot_faces(*mapping_ptr[j], fe_q4, cell, st2); + } + + if (dim>1) + { + ostrstream ost(st2, 99); + ost << "MappingSubface" << dim << "d-" << i << '-' + << mapping_strings[j] << ".dat" << ends; + deallog << st2 << endl; + plot_subfaces(*mapping_ptr[j], fe_q4, cell, st2); + } + + + // Test for transform_*_to_*_cell + if (dim==2 && true) + { + Mapping &mapping=*mapping_ptr[j]; + Point p_unit(6/7.,4/7.); + Point p_real=mapping.transform_unit_to_real_cell(cell, p_unit); + Point p_re_unit=mapping.transform_real_to_unit_cell(cell, p_real); + deallog << "p_unit=" << p_unit << ", p_real=" << p_real + << ", p_re_unit=" << p_re_unit << endl; + } + + delete[] st2; + } + } + + + // delete all triangulations and + // boundary objects + for (unsigned int i=0; i(); + + + // ----------------------- + // Tests for dim=2 + // ----------------------- + mapping_test<2>(); + + + // ----------------------- + // Tests for dim=3 + // ----------------------- + mapping_test<3>(); +} + diff --git a/tests/fe/shapes.cc b/tests/fe/shapes.cc new file mode 100644 index 0000000000..2b79f0925d --- /dev/null +++ b/tests/fe/shapes.cc @@ -0,0 +1,241 @@ +// $Id$ +// (c) Guido Kanschat +// +// Show the shape functions implemented. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +char fname[50]; + +template +inline void +plot_shape_functions(Mapping& mapping, + FiniteElement& finel, + const char* name) +{ + Triangulation tr; + DoFHandler dof(tr); + GridGenerator::hyper_cube(tr, 0., 1.); + DoFHandler::cell_iterator c = dof.begin(); + dof.distribute_dofs(finel); + + const unsigned int div = 11; + + QTrapez<1> q_trapez; + QIterated q(q_trapez, div); + FEValues fe(mapping, finel, q, UpdateFlags(update_values)); + fe.reinit(c); + + sprintf(fname, "Shapes%dd-%s.dat", dim, name); + ofstream gnuplot(fname); + gnuplot.setf(ios::fixed); + gnuplot.precision (2); + + unsigned int k=0; + for (unsigned int mz=0;mz<=((dim>2) ? div : 0) ;++mz) + { + for (unsigned int my=0;my<=((dim>1) ? div : 0) ;++my) + { + for (unsigned int mx=0;mx<=div;++mx) + { + gnuplot << q.point(k); + + for (unsigned int i=0;i +inline void +plot_face_shape_functions(Mapping& mapping, + FiniteElement& finel, + const char* name) +{ + Triangulation tr; + DoFHandler dof(tr); + GridGenerator::hyper_cube(tr, 0., 1.); + tr.refine_global(1); + DoFHandler::active_cell_iterator c = dof.begin_active(); + ++c; + c->set_refine_flag(); + tr.execute_coarsening_and_refinement (); + c = dof.begin_active(); + + dof.distribute_dofs(finel); + + const unsigned int div = 4; + + QTrapez<1> q_trapez; + QIterated q(q_trapez, div); + FEFaceValues fe(mapping, finel, q, UpdateFlags(update_values + | update_q_points)); + FESubfaceValues sub(mapping, finel, q, UpdateFlags(update_values + | update_q_points)); + + sprintf(fname, "ShapesFace%dd-%s.dat", dim, name); + ofstream gnuplot(fname); + gnuplot.setf(ios::fixed); + gnuplot.precision (2); + + for (unsigned int f=0; f::faces_per_cell; ++f) + { + if (!c->face(f)->has_children()) + { + fe.reinit(c, f); + + unsigned int k=0; + for (unsigned int my=0;my<=((dim>2) ? div : 0) ;++my) + { + for (unsigned int mx=0;mx<=div;++mx) + { + gnuplot << fe.quadrature_point(k); + + for (unsigned int i=0;i::subfaces_per_face; ++s) + { + sub.reinit(c, f, s); + + unsigned int k=0; + for (unsigned int my=0;my<=((dim>2) ? div : 0) ;++my) + { + for (unsigned int mx=0;mx<=div;++mx) + { + gnuplot << sub.quadrature_point(k); + + for (unsigned int i=0;i +void plot_face_shape_functions (Mapping<1>&, + FiniteElement<1>&, + const char*) +{} + + +template +void plot_FE_Q_shape_functions() +{ + MappingQ1 m; + FE_Q q1(1); + plot_shape_functions(m, q1, "Q1"); + plot_face_shape_functions(m, q1, "Q1"); + FE_Q q2(2); + plot_shape_functions(m, q2, "Q2"); + plot_face_shape_functions(m, q2, "Q2"); + FE_Q q3(3); + plot_shape_functions(m, q3, "Q3"); + plot_face_shape_functions(m, q3, "Q3"); + FE_Q q4(4); + plot_shape_functions(m, q4, "Q4"); + plot_face_shape_functions(m, q4, "Q4"); +// FE_Q q5(5); +// plot_shape_functions(m, q5, "Q5"); +// FE_Q q6(6); +// plot_shape_functions(m, q6, "Q6"); +// FE_Q q7(7); +// plot_shape_functions(m, q7, "Q7"); +// FE_Q q8(8); +// plot_shape_functions(m, q8, "Q8"); +// FE_Q q9(9); +// plot_shape_functions(m, q9, "Q9"); +// FE_Q q10(10); +// plot_shape_functions(m, q10, "Q10"); +} + + +template +void plot_FE_DGQ_shape_functions() +{ + MappingQ1 m; + FE_DGQ q1(1); + plot_shape_functions(m, q1, "DGQ1"); + plot_face_shape_functions(m, q1, "DGQ1"); + FE_DGQ q2(2); + plot_shape_functions(m, q2, "DGQ2"); + plot_face_shape_functions(m, q2, "DGQ2"); + FE_DGQ q3(3); + plot_shape_functions(m, q3, "DGQ3"); + plot_face_shape_functions(m, q3, "DGQ3"); + FE_DGQ q4(4); + plot_shape_functions(m, q4, "DGQ4"); + plot_face_shape_functions(m, q4, "DGQ4"); +// FE_DGQ q5(5); +// plot_shape_functions(m, q5, "DGQ5"); +// FE_DGQ q6(6); +// plot_shape_functions(m, q6, "DGQ6"); +// FE_DGQ q7(7); +// plot_shape_functions(m, q7, "DGQ7"); +// FE_DGQ q8(8); +// plot_shape_functions(m, q8, "DGQ8"); +// FE_DGQ q9(9); +// plot_shape_functions(m, q9, "DGQ9"); +// FE_DGQ q10(10); +// plot_shape_functions(m, q10, "DGQ10"); +} + + +int +main() +{ + plot_FE_Q_shape_functions<1>(); + plot_FE_Q_shape_functions<2>(); + plot_FE_DGQ_shape_functions<2>(); +// plot_FE_Q_shape_functions<3>(); + + // FESystem test. + MappingQ1<2> m; + FESystem<2> q2_q3(FE_Q<2>(2), 1, + FE_Q<2>(3), 1); +// plot_shape_functions(m, q2_q3, "Q2_Q3"); + + return 0; +} + + + -- 2.39.5