From: Wolfgang Bangerth Date: Thu, 2 Apr 1998 15:05:33 +0000 (+0000) Subject: Rewrite of the poisson example, using the ParameterHandler class. X-Git-Tag: v8.0.0~23132 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0f9b7a90c683421ff70000dbbc05e2d55f9c4618;p=dealii.git Rewrite of the poisson example, using the ParameterHandler class. git-svn-id: https://svn.dealii.org/trunk@122 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/Attic/examples/poisson/Makefile b/deal.II/deal.II/Attic/examples/poisson/Makefile new file mode 100644 index 0000000000..91ad7869d0 --- /dev/null +++ b/deal.II/deal.II/Attic/examples/poisson/Makefile @@ -0,0 +1,59 @@ +# $Id$ + +CXX = c++ +INCLUDE = -I../../include -I../../../lac/include \ + -I../../../base/include -I../../../mia/include +CXXFLAGS.g= -DDEBUG -g -Wall -pedantic -Wconversion \ + -Winline -Woverloaded-virtual $(INCLUDE) +CXXFLAGS =-O3 -Wuninitialized -finline-functions -ffast-math \ + -funroll-loops -felide-constructors -fnonnull-objects $(INCLUDE) + +ifeq ($(shell uname),SunOS) +#CXXFLAGS := -V2.7.2.3 $(CXXFLAGS) +#CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g) +INCLUDE2.8 = -I/usr/local/source/libstdc++-2.8.0/libstdc++ \ + -I/usr/local/source/libstdc++-2.8.0/libstdc++/stl \ + -I/usr/local/source/libstdc++-2.8.0/libio +CXXFLAGS := $(CXXFLAGS) -fno-rtti -fno-exceptions $(INCLUDE2.8) +CXXFLAGS.g:= $(CXXFLAGS.g) -W -fno-rtti -fno-exceptions $(INCLUDE2.8) +endif + +LIBPATH = -L../../lib -L../../../lac/lib +LIBFILES.g= ../../lib/libbasic.g.a ../../lib/libgrid.g.a \ + ../../lib/libfe.g.a ../../lib/libnumerics.g.a \ + ../../../lac/lib/liblac.g.a +LIBFILES = $(LIBFILES.g:.g.a=.a) +LIBS.g = $(LIBPATH) -lbasic.g -lgrid.g -lfe.g -lnumerics.g -llac.g +LIBS = $(LIBS.g:.g=) + + +cc-files = $(wildcard *.cc) +o-files = $(cc-files:.cc=.o) +go-files = $(cc-files:.cc=.go) +h-files = $(wildcard ../include/*/*.h) + + + + +%.go : %.cc poisson.h + @echo ============================ Compiling with debugging information: $< + @echo $(CXX) ... -c $< -o $@ + @$(CXX) $(CXXFLAGS.g) -c $< -o $@ +%.o : %.cc poisson.h + @echo ============================ Compiling with optimization: $< + @echo $(CXX) ... -c $< -o $@ + @$(CXX) $(CXXFLAGS) -c $< -o $@ + + +poisson : $(go-files) $(LIBFILES.g) + @echo ============================ Linking poisson + $(CXX) -g $(CXXFLAGS.g) $(go-files) -o $@ ../../../mia/control.o $(LIBS.g) -lg++ + + + +clean: + rm -f *.o *.go *~ poisson + + + +.PHONY: clean diff --git a/deal.II/deal.II/Attic/examples/poisson/equation.cc b/deal.II/deal.II/Attic/examples/poisson/equation.cc new file mode 100644 index 0000000000..c995c110e0 --- /dev/null +++ b/deal.II/deal.II/Attic/examples/poisson/equation.cc @@ -0,0 +1,91 @@ +/* $Id$ */ + +#include "poisson.h" + + + +template +double RightHandSide::operator () (const Point &p) const { + const double pi = 3.1415926536; + switch (dim) + { + case 1: + return p(0)*p(0)*p(0)-3./2.*p(0)*p(0)-6*p(0)+3; + case 2: + return (-2.0*cos(pi*p(0)/2)*p(1)*sin(pi*p(1)) + + 2.0*p(0)*sin(pi*p(0)/2)*pi*p(1)*sin(pi*p(1)) + + 5.0/4.0*p(0)*p(0)*cos(pi*p(0)/2)*pi*pi*p(1)*sin(pi*p(1)) - + 2.0*p(0)*p(0)*cos(pi*p(0)/2)*cos(pi*p(1))*pi); + default: + return 0; + }; +}; + + + + + +void PoissonEquation<1>::assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues<1> &fe_values, + const Triangulation<1>::cell_iterator &) const { + for (unsigned int point=0; point::assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues<2> &fe_values, + const Triangulation<2>::cell_iterator &) const { + for (unsigned int point=0; point +void PoissonEquation::assemble (dFMatrix &, + const FEValues &, + const Triangulation::cell_iterator &) const { + Assert (false, ExcPureVirtualFunctionCalled()); +}; + + + +template +void PoissonEquation::assemble (dVector &, + const FEValues &, + const Triangulation::cell_iterator &) const { + Assert (false, ExcPureVirtualFunctionCalled()); +}; + + + + + +template class RightHandSide<1>; +template class RightHandSide<2>; + + +template class PoissonEquation<1>; +template class PoissonEquation<2>; diff --git a/deal.II/deal.II/Attic/examples/poisson/poisson.cc b/deal.II/deal.II/Attic/examples/poisson/poisson.cc index a9987b15b3..c2fb186f3f 100644 --- a/deal.II/deal.II/Attic/examples/poisson/poisson.cc +++ b/deal.II/deal.II/Attic/examples/poisson/poisson.cc @@ -1,196 +1,23 @@ /* $Id$ */ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include "poisson.h" -#include -#include -#include -extern "C" { -# include -} -extern TriaActiveIterator<1,CellAccessor<1> > x; -extern TriaActiveIterator<2,CellAccessor<2> > y; -extern TriaRawIterator<1,DoFLineAccessor<1,LineAccessor<1> > > z; - - - -template -class PoissonEquation : public Equation { - public: - PoissonEquation () : - Equation(1) {}; - - virtual void assemble (dFMatrix &cell_matrix, - dVector &rhs, - const FEValues &fe_values, - const Triangulation::cell_iterator &cell) const; - double right_hand_side (const Point &) const; -}; - - - -template -inline -double PoissonEquation::right_hand_side (const Point &p) const { - const double pi = 3.1415926536; - switch (dim) - { - case 1: -// return ((1-4*3.1415926536*3.1415926536) * -// cos(2*3.1415926536*p(0))); - return p(0)*p(0)*p(0)-3./2.*p(0)*p(0)-6*p(0)+3; - case 2: -// return ((1-3.1415926536*3.1415926536) * -// cos(3.1415926536*p(0)) * -// cos(3.1415926536*p(1))); -// return (p(0)*p(0)*p(0)+p(1)*p(1)*p(1) -// - 3./2.*(p(0)*p(0)+p(1)*p(1)) -// - 6*(p(0)+p(1)) -// + 6); - return (-2.0*cos(pi*p(0)/2)*p(1)*sin(pi*p(1)) + - 2.0*p(0)*sin(pi*p(0)/2)*pi*p(1)*sin(pi*p(1)) + - 5.0/4.0*p(0)*p(0)*cos(pi*p(0)/2)*pi*pi*p(1)*sin(pi*p(1)) - - 2.0*p(0)*p(0)*cos(pi*p(0)/2)*cos(pi*p(1))*pi); - default: - return 0; - }; -}; - - - -void PoissonEquation<1>::assemble (dFMatrix &cell_matrix, - dVector &rhs, - const FEValues<1> &fe_values, - const Triangulation<1>::cell_iterator &) const { - for (unsigned int point=0; point::assemble (dFMatrix &cell_matrix, - dVector &rhs, - const FEValues<2> &fe_values, - const Triangulation<2>::cell_iterator &) const { - for (unsigned int point=0; point tria; - DoFHandler<2> dof(&tria); - FELinear<2> fe; - ProblemBase<2> problem; - problem.set_tria_and_dof (&tria, &dof); - PoissonEquation<2> equation; - QGauss4<2> quadrature; - - - cout << "Making grid..." << endl; - - tria.create_hypercube (); -// HyperBallBoundary<2> boundary(Point<2>(2,3), 4); -// tria.create_hyper_ball(Point<2>(2,3),4); -// tria.set_boundary (&boundary); - -// tria.refine_global (1); -// (--tria.last_active())->set_refine_flag(); -// tria.execute_refinement (); -// tria.begin_active(2)->set_refine_flag(); -// tria.execute_refinement (); -// tria.refine_global (5); - - - const unsigned int dim=2; - tria.refine_global (1); - - Triangulation::active_cell_iterator cell, endc; - for (int i=0; i<12; ++i) +int main (int argc, char **argv) { + if (argc!=2) { - int n_levels = tria.n_levels(); - cell = tria.begin_active(); - endc = tria.end(); - - for (; cell!=endc; ++cell) - { - double r = rand()*1.0/RAND_MAX, - weight = 1.* - (cell->level()*cell->level()) / - (n_levels*n_levels); - - if (r <= 0.5*weight) - cell->set_refine_flag (); - }; - - tria.execute_refinement (); + cout << "Usage: poisson parameterfile" << endl << endl; + return 1; }; - tria.refine_global (1); - - - cout << "Distributing dofs... "; - dof.distribute_dofs (fe); - cout << dof.n_dofs() << " degrees of freedom." << endl; - - cout << "Assembling matrices..." << endl; - FEValues<2>::UpdateStruct update_flags; - update_flags.q_points = update_flags.gradients = true; - update_flags.jacobians = update_flags.JxW_values = true; - - ProblemBase<2>::DirichletBC dirichlet_bc; - ZeroFunction<2> zero; - dirichlet_bc[0] = &zero; - problem.assemble (equation, quadrature, fe, update_flags, dirichlet_bc); - cout << "Solving..." << endl; - problem.solve (); + PoissonProblem<2> poisson; + MultipleParameterLoop input_data; - cout << "Printing..." << endl; - DataOut<2> out; - ofstream gnuplot("gnuplot.out.5"); - problem.fill_data (out); - out.write_gnuplot (gnuplot); - gnuplot.close (); + poisson.declare_parameters(input_data); + input_data.read_input (argv[1]); + input_data.loop (poisson); return 0; }; diff --git a/deal.II/deal.II/Attic/examples/poisson/poisson.h b/deal.II/deal.II/Attic/examples/poisson/poisson.h new file mode 100644 index 0000000000..d41e6bf432 --- /dev/null +++ b/deal.II/deal.II/Attic/examples/poisson/poisson.h @@ -0,0 +1,109 @@ +/*---------------------------- poisson.h ---------------------------*/ +/* $Id$ */ +#ifndef __poisson_H +#define __poisson_H +/*---------------------------- poisson.h ---------------------------*/ + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include +#include +#include +extern "C" { +# include +} + + + + + + +template +class RightHandSide : public Function { + public: + /** + * Return the value of the function + * at the given point. + */ + virtual double operator () (const Point &p) const; +}; + + + + + + +template +class PoissonEquation : public Equation { + public: + PoissonEquation (const Function &rhs) : + Equation(1), + right_hand_side (rhs) {}; + + virtual void assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + virtual void assemble (dFMatrix &cell_matrix, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + virtual void assemble (dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + protected: + const Function &right_hand_side; +}; + + + + + + +template +class PoissonProblem : public ProblemBase, + public MultipleParameterLoop::UserClass { + public: + PoissonProblem (); + + void clear (); + + virtual void create_new (const unsigned int run_no); + virtual void declare_parameters (ParameterHandler &prm); + virtual void run (ParameterHandler &prm); + + + bool make_grid (ParameterHandler &prm); + void make_ball_grid (); + void make_curved_line_grid (); + void make_zoom_in_grid (); + void make_random_grid (); + + protected: + Triangulation *tria; + DoFHandler *dof; +}; + + + + + +/*---------------------------- poisson.h ---------------------------*/ +/* end of #ifndef __poisson_H */ +#endif +/*---------------------------- poisson.h ---------------------------*/ diff --git a/deal.II/deal.II/Attic/examples/poisson/problem.cc b/deal.II/deal.II/Attic/examples/poisson/problem.cc new file mode 100644 index 0000000000..acef23aa40 --- /dev/null +++ b/deal.II/deal.II/Attic/examples/poisson/problem.cc @@ -0,0 +1,248 @@ +/* $Id$ */ + + +#include "poisson.h" + + + + +template +class CurvedLine : + public StraightBoundary { + public: + virtual Point in_between (const PointArray &neighbors) const { + Point middle = StraightBoundary::in_between(neighbors); + double x=middle(0), + y=middle(1); + + if (y +PoissonProblem::PoissonProblem () : + tria(0), dof(0) {}; + + + + +template +void PoissonProblem::clear () { + if (tria != 0) delete tria; + if (dof != 0) delete dof; + + ProblemBase::clear (); +}; + + + + +template +void PoissonProblem::create_new (const unsigned int) { + clear (); + + tria = new Triangulation(); + dof = new DoFHandler (tria); + set_tria_and_dof (tria, dof); +}; + + + + +template +void PoissonProblem::declare_parameters (ParameterHandler &prm) { + if (dim>=2) + prm.declare_entry ("Test run", "zoom in", "zoom in\\|ball\\|curved line\\|random"); + else + prm.declare_entry ("Test run", "zoom in", "zoom in\\|random"); + + prm.declare_entry ("Global refinement", "0", ParameterHandler::RegularExpressions::Integer); +}; + + + + +template +bool PoissonProblem::make_grid (ParameterHandler &prm) { + String test = prm.get ("Test run"); + unsigned int test_case; + if (test=="zoom in") test_case = 1; + else + if (test=="ball") test_case = 2; + else + if (test=="curved line") test_case = 3; + else + if (test=="random") test_case = 4; + else + cerr << "This test seems not to be implemented!" << endl; + + switch (test_case) + { + case 1: + make_zoom_in_grid (); + return true; + case 2: + make_ball_grid (); + return true; + case 3: + make_curved_line_grid (); + return true; + case 4: + make_random_grid (); + return true; + default: + return false; + }; + + int refine_global = prm.get_integer ("Global refinement"); + if ((refine_global < 0) || (refine_global>10)) + return false; + else + tria->refine_global (refine_global); +}; + + + + +template +void PoissonProblem::make_zoom_in_grid () { + tria->create_hypercube (); + // refine first cell + tria->begin_active()->set_refine_flag(); + tria->execute_refinement (); + // refine first active cell + // on coarsest level + tria->begin_active()->set_refine_flag (); + tria->execute_refinement (); + + Triangulation::active_cell_iterator cell; + for (int i=0; i<17; ++i) + { + // refine the presently + // second last cell 17 + // times + cell = tria->last_active(tria->n_levels()-1); + --cell; + cell->set_refine_flag (); + tria->execute_refinement (); + }; +}; + + + +template +void PoissonProblem::make_ball_grid () { + // make ball grid around origin with + // unit radius + const Point origin; + static const HyperBallBoundary boundary(origin, 1.); + + tria->create_hyper_ball (origin, 1.); + tria->set_boundary (&boundary); +}; + + + + +template +void PoissonProblem::make_curved_line_grid () { + // set the boundary function + static const CurvedLine boundary; + + tria->create_hypercube (); + tria->set_boundary (&boundary); +}; + + + +template +void PoissonProblem::make_random_grid () { + tria->refine_global (1); + + Triangulation::active_cell_iterator cell, endc; + for (int i=0; i<12; ++i) + { + int n_levels = tria->n_levels(); + cell = tria->begin_active(); + endc = tria->end(); + + for (; cell!=endc; ++cell) + { + double r = rand()*1.0/RAND_MAX, + weight = 1.* + (cell->level()*cell->level()) / + (n_levels*n_levels); + + if (r <= 0.5*weight) + cell->set_refine_flag (); + }; + + tria->execute_refinement (); + }; +}; + + + + + +template +void PoissonProblem::run (ParameterHandler &prm) { + cout << "Test case = " << prm.get ("Test run") << endl + << endl; + + cout << " Making grid..." << endl; + if (!make_grid (prm)) + return; + + + FELinear fe; + static const RightHandSide rhs; + PoissonEquation equation (rhs); + QGauss4 quadrature; + + cout << "Distributing dofs... "; + dof->distribute_dofs (fe); + cout << dof->n_dofs() << " degrees of freedom." << endl; + + cout << "Assembling matrices..." << endl; + FEValues::UpdateStruct update_flags; + update_flags.q_points = update_flags.gradients = true; + update_flags.jacobians = update_flags.JxW_values = true; + + ProblemBase::DirichletBC dirichlet_bc; + ZeroFunction zero; + dirichlet_bc[0] = &zero; + assemble (equation, quadrature, fe, update_flags, dirichlet_bc); + + cout << "Solving..." << endl; + solve (); + + cout << "Printing..." << endl; + DataOut out; + ofstream gnuplot("gnuplot.out.5"); + fill_data (out); + out.write_gnuplot (gnuplot); + gnuplot.close (); +}; + + + + + +template class PoissonProblem<1>; +template class PoissonProblem<2>; diff --git a/tests/big-tests/poisson/Makefile b/tests/big-tests/poisson/Makefile new file mode 100644 index 0000000000..91ad7869d0 --- /dev/null +++ b/tests/big-tests/poisson/Makefile @@ -0,0 +1,59 @@ +# $Id$ + +CXX = c++ +INCLUDE = -I../../include -I../../../lac/include \ + -I../../../base/include -I../../../mia/include +CXXFLAGS.g= -DDEBUG -g -Wall -pedantic -Wconversion \ + -Winline -Woverloaded-virtual $(INCLUDE) +CXXFLAGS =-O3 -Wuninitialized -finline-functions -ffast-math \ + -funroll-loops -felide-constructors -fnonnull-objects $(INCLUDE) + +ifeq ($(shell uname),SunOS) +#CXXFLAGS := -V2.7.2.3 $(CXXFLAGS) +#CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g) +INCLUDE2.8 = -I/usr/local/source/libstdc++-2.8.0/libstdc++ \ + -I/usr/local/source/libstdc++-2.8.0/libstdc++/stl \ + -I/usr/local/source/libstdc++-2.8.0/libio +CXXFLAGS := $(CXXFLAGS) -fno-rtti -fno-exceptions $(INCLUDE2.8) +CXXFLAGS.g:= $(CXXFLAGS.g) -W -fno-rtti -fno-exceptions $(INCLUDE2.8) +endif + +LIBPATH = -L../../lib -L../../../lac/lib +LIBFILES.g= ../../lib/libbasic.g.a ../../lib/libgrid.g.a \ + ../../lib/libfe.g.a ../../lib/libnumerics.g.a \ + ../../../lac/lib/liblac.g.a +LIBFILES = $(LIBFILES.g:.g.a=.a) +LIBS.g = $(LIBPATH) -lbasic.g -lgrid.g -lfe.g -lnumerics.g -llac.g +LIBS = $(LIBS.g:.g=) + + +cc-files = $(wildcard *.cc) +o-files = $(cc-files:.cc=.o) +go-files = $(cc-files:.cc=.go) +h-files = $(wildcard ../include/*/*.h) + + + + +%.go : %.cc poisson.h + @echo ============================ Compiling with debugging information: $< + @echo $(CXX) ... -c $< -o $@ + @$(CXX) $(CXXFLAGS.g) -c $< -o $@ +%.o : %.cc poisson.h + @echo ============================ Compiling with optimization: $< + @echo $(CXX) ... -c $< -o $@ + @$(CXX) $(CXXFLAGS) -c $< -o $@ + + +poisson : $(go-files) $(LIBFILES.g) + @echo ============================ Linking poisson + $(CXX) -g $(CXXFLAGS.g) $(go-files) -o $@ ../../../mia/control.o $(LIBS.g) -lg++ + + + +clean: + rm -f *.o *.go *~ poisson + + + +.PHONY: clean diff --git a/tests/big-tests/poisson/equation.cc b/tests/big-tests/poisson/equation.cc new file mode 100644 index 0000000000..c995c110e0 --- /dev/null +++ b/tests/big-tests/poisson/equation.cc @@ -0,0 +1,91 @@ +/* $Id$ */ + +#include "poisson.h" + + + +template +double RightHandSide::operator () (const Point &p) const { + const double pi = 3.1415926536; + switch (dim) + { + case 1: + return p(0)*p(0)*p(0)-3./2.*p(0)*p(0)-6*p(0)+3; + case 2: + return (-2.0*cos(pi*p(0)/2)*p(1)*sin(pi*p(1)) + + 2.0*p(0)*sin(pi*p(0)/2)*pi*p(1)*sin(pi*p(1)) + + 5.0/4.0*p(0)*p(0)*cos(pi*p(0)/2)*pi*pi*p(1)*sin(pi*p(1)) - + 2.0*p(0)*p(0)*cos(pi*p(0)/2)*cos(pi*p(1))*pi); + default: + return 0; + }; +}; + + + + + +void PoissonEquation<1>::assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues<1> &fe_values, + const Triangulation<1>::cell_iterator &) const { + for (unsigned int point=0; point::assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues<2> &fe_values, + const Triangulation<2>::cell_iterator &) const { + for (unsigned int point=0; point +void PoissonEquation::assemble (dFMatrix &, + const FEValues &, + const Triangulation::cell_iterator &) const { + Assert (false, ExcPureVirtualFunctionCalled()); +}; + + + +template +void PoissonEquation::assemble (dVector &, + const FEValues &, + const Triangulation::cell_iterator &) const { + Assert (false, ExcPureVirtualFunctionCalled()); +}; + + + + + +template class RightHandSide<1>; +template class RightHandSide<2>; + + +template class PoissonEquation<1>; +template class PoissonEquation<2>; diff --git a/tests/big-tests/poisson/poisson.cc b/tests/big-tests/poisson/poisson.cc index a9987b15b3..c2fb186f3f 100644 --- a/tests/big-tests/poisson/poisson.cc +++ b/tests/big-tests/poisson/poisson.cc @@ -1,196 +1,23 @@ /* $Id$ */ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include "poisson.h" -#include -#include -#include -extern "C" { -# include -} -extern TriaActiveIterator<1,CellAccessor<1> > x; -extern TriaActiveIterator<2,CellAccessor<2> > y; -extern TriaRawIterator<1,DoFLineAccessor<1,LineAccessor<1> > > z; - - - -template -class PoissonEquation : public Equation { - public: - PoissonEquation () : - Equation(1) {}; - - virtual void assemble (dFMatrix &cell_matrix, - dVector &rhs, - const FEValues &fe_values, - const Triangulation::cell_iterator &cell) const; - double right_hand_side (const Point &) const; -}; - - - -template -inline -double PoissonEquation::right_hand_side (const Point &p) const { - const double pi = 3.1415926536; - switch (dim) - { - case 1: -// return ((1-4*3.1415926536*3.1415926536) * -// cos(2*3.1415926536*p(0))); - return p(0)*p(0)*p(0)-3./2.*p(0)*p(0)-6*p(0)+3; - case 2: -// return ((1-3.1415926536*3.1415926536) * -// cos(3.1415926536*p(0)) * -// cos(3.1415926536*p(1))); -// return (p(0)*p(0)*p(0)+p(1)*p(1)*p(1) -// - 3./2.*(p(0)*p(0)+p(1)*p(1)) -// - 6*(p(0)+p(1)) -// + 6); - return (-2.0*cos(pi*p(0)/2)*p(1)*sin(pi*p(1)) + - 2.0*p(0)*sin(pi*p(0)/2)*pi*p(1)*sin(pi*p(1)) + - 5.0/4.0*p(0)*p(0)*cos(pi*p(0)/2)*pi*pi*p(1)*sin(pi*p(1)) - - 2.0*p(0)*p(0)*cos(pi*p(0)/2)*cos(pi*p(1))*pi); - default: - return 0; - }; -}; - - - -void PoissonEquation<1>::assemble (dFMatrix &cell_matrix, - dVector &rhs, - const FEValues<1> &fe_values, - const Triangulation<1>::cell_iterator &) const { - for (unsigned int point=0; point::assemble (dFMatrix &cell_matrix, - dVector &rhs, - const FEValues<2> &fe_values, - const Triangulation<2>::cell_iterator &) const { - for (unsigned int point=0; point tria; - DoFHandler<2> dof(&tria); - FELinear<2> fe; - ProblemBase<2> problem; - problem.set_tria_and_dof (&tria, &dof); - PoissonEquation<2> equation; - QGauss4<2> quadrature; - - - cout << "Making grid..." << endl; - - tria.create_hypercube (); -// HyperBallBoundary<2> boundary(Point<2>(2,3), 4); -// tria.create_hyper_ball(Point<2>(2,3),4); -// tria.set_boundary (&boundary); - -// tria.refine_global (1); -// (--tria.last_active())->set_refine_flag(); -// tria.execute_refinement (); -// tria.begin_active(2)->set_refine_flag(); -// tria.execute_refinement (); -// tria.refine_global (5); - - - const unsigned int dim=2; - tria.refine_global (1); - - Triangulation::active_cell_iterator cell, endc; - for (int i=0; i<12; ++i) +int main (int argc, char **argv) { + if (argc!=2) { - int n_levels = tria.n_levels(); - cell = tria.begin_active(); - endc = tria.end(); - - for (; cell!=endc; ++cell) - { - double r = rand()*1.0/RAND_MAX, - weight = 1.* - (cell->level()*cell->level()) / - (n_levels*n_levels); - - if (r <= 0.5*weight) - cell->set_refine_flag (); - }; - - tria.execute_refinement (); + cout << "Usage: poisson parameterfile" << endl << endl; + return 1; }; - tria.refine_global (1); - - - cout << "Distributing dofs... "; - dof.distribute_dofs (fe); - cout << dof.n_dofs() << " degrees of freedom." << endl; - - cout << "Assembling matrices..." << endl; - FEValues<2>::UpdateStruct update_flags; - update_flags.q_points = update_flags.gradients = true; - update_flags.jacobians = update_flags.JxW_values = true; - - ProblemBase<2>::DirichletBC dirichlet_bc; - ZeroFunction<2> zero; - dirichlet_bc[0] = &zero; - problem.assemble (equation, quadrature, fe, update_flags, dirichlet_bc); - cout << "Solving..." << endl; - problem.solve (); + PoissonProblem<2> poisson; + MultipleParameterLoop input_data; - cout << "Printing..." << endl; - DataOut<2> out; - ofstream gnuplot("gnuplot.out.5"); - problem.fill_data (out); - out.write_gnuplot (gnuplot); - gnuplot.close (); + poisson.declare_parameters(input_data); + input_data.read_input (argv[1]); + input_data.loop (poisson); return 0; }; diff --git a/tests/big-tests/poisson/poisson.h b/tests/big-tests/poisson/poisson.h new file mode 100644 index 0000000000..d41e6bf432 --- /dev/null +++ b/tests/big-tests/poisson/poisson.h @@ -0,0 +1,109 @@ +/*---------------------------- poisson.h ---------------------------*/ +/* $Id$ */ +#ifndef __poisson_H +#define __poisson_H +/*---------------------------- poisson.h ---------------------------*/ + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include +#include +#include +extern "C" { +# include +} + + + + + + +template +class RightHandSide : public Function { + public: + /** + * Return the value of the function + * at the given point. + */ + virtual double operator () (const Point &p) const; +}; + + + + + + +template +class PoissonEquation : public Equation { + public: + PoissonEquation (const Function &rhs) : + Equation(1), + right_hand_side (rhs) {}; + + virtual void assemble (dFMatrix &cell_matrix, + dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + virtual void assemble (dFMatrix &cell_matrix, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + virtual void assemble (dVector &rhs, + const FEValues &fe_values, + const Triangulation::cell_iterator &cell) const; + protected: + const Function &right_hand_side; +}; + + + + + + +template +class PoissonProblem : public ProblemBase, + public MultipleParameterLoop::UserClass { + public: + PoissonProblem (); + + void clear (); + + virtual void create_new (const unsigned int run_no); + virtual void declare_parameters (ParameterHandler &prm); + virtual void run (ParameterHandler &prm); + + + bool make_grid (ParameterHandler &prm); + void make_ball_grid (); + void make_curved_line_grid (); + void make_zoom_in_grid (); + void make_random_grid (); + + protected: + Triangulation *tria; + DoFHandler *dof; +}; + + + + + +/*---------------------------- poisson.h ---------------------------*/ +/* end of #ifndef __poisson_H */ +#endif +/*---------------------------- poisson.h ---------------------------*/ diff --git a/tests/big-tests/poisson/problem.cc b/tests/big-tests/poisson/problem.cc new file mode 100644 index 0000000000..acef23aa40 --- /dev/null +++ b/tests/big-tests/poisson/problem.cc @@ -0,0 +1,248 @@ +/* $Id$ */ + + +#include "poisson.h" + + + + +template +class CurvedLine : + public StraightBoundary { + public: + virtual Point in_between (const PointArray &neighbors) const { + Point middle = StraightBoundary::in_between(neighbors); + double x=middle(0), + y=middle(1); + + if (y +PoissonProblem::PoissonProblem () : + tria(0), dof(0) {}; + + + + +template +void PoissonProblem::clear () { + if (tria != 0) delete tria; + if (dof != 0) delete dof; + + ProblemBase::clear (); +}; + + + + +template +void PoissonProblem::create_new (const unsigned int) { + clear (); + + tria = new Triangulation(); + dof = new DoFHandler (tria); + set_tria_and_dof (tria, dof); +}; + + + + +template +void PoissonProblem::declare_parameters (ParameterHandler &prm) { + if (dim>=2) + prm.declare_entry ("Test run", "zoom in", "zoom in\\|ball\\|curved line\\|random"); + else + prm.declare_entry ("Test run", "zoom in", "zoom in\\|random"); + + prm.declare_entry ("Global refinement", "0", ParameterHandler::RegularExpressions::Integer); +}; + + + + +template +bool PoissonProblem::make_grid (ParameterHandler &prm) { + String test = prm.get ("Test run"); + unsigned int test_case; + if (test=="zoom in") test_case = 1; + else + if (test=="ball") test_case = 2; + else + if (test=="curved line") test_case = 3; + else + if (test=="random") test_case = 4; + else + cerr << "This test seems not to be implemented!" << endl; + + switch (test_case) + { + case 1: + make_zoom_in_grid (); + return true; + case 2: + make_ball_grid (); + return true; + case 3: + make_curved_line_grid (); + return true; + case 4: + make_random_grid (); + return true; + default: + return false; + }; + + int refine_global = prm.get_integer ("Global refinement"); + if ((refine_global < 0) || (refine_global>10)) + return false; + else + tria->refine_global (refine_global); +}; + + + + +template +void PoissonProblem::make_zoom_in_grid () { + tria->create_hypercube (); + // refine first cell + tria->begin_active()->set_refine_flag(); + tria->execute_refinement (); + // refine first active cell + // on coarsest level + tria->begin_active()->set_refine_flag (); + tria->execute_refinement (); + + Triangulation::active_cell_iterator cell; + for (int i=0; i<17; ++i) + { + // refine the presently + // second last cell 17 + // times + cell = tria->last_active(tria->n_levels()-1); + --cell; + cell->set_refine_flag (); + tria->execute_refinement (); + }; +}; + + + +template +void PoissonProblem::make_ball_grid () { + // make ball grid around origin with + // unit radius + const Point origin; + static const HyperBallBoundary boundary(origin, 1.); + + tria->create_hyper_ball (origin, 1.); + tria->set_boundary (&boundary); +}; + + + + +template +void PoissonProblem::make_curved_line_grid () { + // set the boundary function + static const CurvedLine boundary; + + tria->create_hypercube (); + tria->set_boundary (&boundary); +}; + + + +template +void PoissonProblem::make_random_grid () { + tria->refine_global (1); + + Triangulation::active_cell_iterator cell, endc; + for (int i=0; i<12; ++i) + { + int n_levels = tria->n_levels(); + cell = tria->begin_active(); + endc = tria->end(); + + for (; cell!=endc; ++cell) + { + double r = rand()*1.0/RAND_MAX, + weight = 1.* + (cell->level()*cell->level()) / + (n_levels*n_levels); + + if (r <= 0.5*weight) + cell->set_refine_flag (); + }; + + tria->execute_refinement (); + }; +}; + + + + + +template +void PoissonProblem::run (ParameterHandler &prm) { + cout << "Test case = " << prm.get ("Test run") << endl + << endl; + + cout << " Making grid..." << endl; + if (!make_grid (prm)) + return; + + + FELinear fe; + static const RightHandSide rhs; + PoissonEquation equation (rhs); + QGauss4 quadrature; + + cout << "Distributing dofs... "; + dof->distribute_dofs (fe); + cout << dof->n_dofs() << " degrees of freedom." << endl; + + cout << "Assembling matrices..." << endl; + FEValues::UpdateStruct update_flags; + update_flags.q_points = update_flags.gradients = true; + update_flags.jacobians = update_flags.JxW_values = true; + + ProblemBase::DirichletBC dirichlet_bc; + ZeroFunction zero; + dirichlet_bc[0] = &zero; + assemble (equation, quadrature, fe, update_flags, dirichlet_bc); + + cout << "Solving..." << endl; + solve (); + + cout << "Printing..." << endl; + DataOut out; + ofstream gnuplot("gnuplot.out.5"); + fill_data (out); + out.write_gnuplot (gnuplot); + gnuplot.close (); +}; + + + + + +template class PoissonProblem<1>; +template class PoissonProblem<2>;