--- /dev/null
+# $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
--- /dev/null
+/* $Id$ */
+
+#include "poisson.h"
+
+
+
+template <int dim>
+double RightHandSide<dim>::operator () (const Point<dim> &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<fe_values.n_quadrature_points; ++point)
+ for (unsigned int i=0; i<fe_values.total_dofs; ++i)
+ {
+ for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+ cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
+ fe_values.shape_grad(j,point)) *
+ fe_values.JxW(point);
+ rhs(i) += fe_values.shape_value(i,point) *
+ right_hand_side(fe_values.quadrature_point(point)) *
+ fe_values.JxW(point);
+ };
+};
+
+
+
+void PoissonEquation<2>::assemble (dFMatrix &cell_matrix,
+ dVector &rhs,
+ const FEValues<2> &fe_values,
+ const Triangulation<2>::cell_iterator &) const {
+ for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+ for (unsigned int i=0; i<fe_values.total_dofs; ++i)
+ {
+ for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+ cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
+ fe_values.shape_grad(j,point)) *
+ fe_values.JxW(point);
+ rhs(i) += fe_values.shape_value(i,point) *
+ right_hand_side(fe_values.quadrature_point(point)) *
+ fe_values.JxW(point);
+ };
+};
+
+
+
+template <int dim>
+void PoissonEquation<dim>::assemble (dFMatrix &,
+ const FEValues<dim> &,
+ const Triangulation<dim>::cell_iterator &) const {
+ Assert (false, ExcPureVirtualFunctionCalled());
+};
+
+
+
+template <int dim>
+void PoissonEquation<dim>::assemble (dVector &,
+ const FEValues<dim> &,
+ const Triangulation<dim>::cell_iterator &) const {
+ Assert (false, ExcPureVirtualFunctionCalled());
+};
+
+
+
+
+
+template class RightHandSide<1>;
+template class RightHandSide<2>;
+
+
+template class PoissonEquation<1>;
+template class PoissonEquation<2>;
/* $Id$ */
-#include <grid/tria.h>
-#include <grid/dof.h>
-#include <grid/tria_accessor.h>
-#include <grid/dof_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_boundary.h>
-#include <grid/dof_constraints.h>
-#include <basic/data_io.h>
-#include <basic/function.h>
-#include <fe/fe_lib.h>
-#include <fe/quadrature_lib.h>
-#include <numerics/base.h>
-#include <numerics/assembler.h>
-#include <lac/dsmatrix.h>
+#include "poisson.h"
-#include <map.h>
-#include <fstream.h>
-#include <cmath>
-extern "C" {
-# include <stdlib.h>
-}
-extern TriaActiveIterator<1,CellAccessor<1> > x;
-extern TriaActiveIterator<2,CellAccessor<2> > y;
-extern TriaRawIterator<1,DoFLineAccessor<1,LineAccessor<1> > > z;
-
-
-
-template <int dim>
-class PoissonEquation : public Equation<dim> {
- public:
- PoissonEquation () :
- Equation<dim>(1) {};
-
- virtual void assemble (dFMatrix &cell_matrix,
- dVector &rhs,
- const FEValues<dim> &fe_values,
- const Triangulation<dim>::cell_iterator &cell) const;
- double right_hand_side (const Point<dim> &) const;
-};
-
-
-
-template <int dim>
-inline
-double PoissonEquation<dim>::right_hand_side (const Point<dim> &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<fe_values.n_quadrature_points; ++point)
- for (unsigned int i=0; i<fe_values.total_dofs; ++i)
- {
- for (unsigned int j=0; j<fe_values.total_dofs; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point) +
- fe_values.shape_value(i,point) *
- fe_values.shape_value(j,point)) *
- fe_values.JxW(point);
- rhs(i) += fe_values.shape_value(i,point) *
- right_hand_side(fe_values.quadrature_point(point)) *
- fe_values.JxW(point);
- };
-};
-
-
-
-void PoissonEquation<2>::assemble (dFMatrix &cell_matrix,
- dVector &rhs,
- const FEValues<2> &fe_values,
- const Triangulation<2>::cell_iterator &) const {
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
- for (unsigned int i=0; i<fe_values.total_dofs; ++i)
- {
- for (unsigned int j=0; j<fe_values.total_dofs; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point)/* +
- fe_values.shape_value(i,point) *
- fe_values.shape_value(j,point)*/) *
- fe_values.JxW(point);
- rhs(i) += fe_values.shape_value(i,point) *
- right_hand_side(fe_values.quadrature_point(point)) *
- fe_values.JxW(point);
- };
-};
-
-
-
-
-int main () {
- Triangulation<2> 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<dim>::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;
};
--- /dev/null
+/*---------------------------- poisson.h ---------------------------*/
+/* $Id$ */
+#ifndef __poisson_H
+#define __poisson_H
+/*---------------------------- poisson.h ---------------------------*/
+
+
+#include <grid/tria.h>
+#include <grid/dof.h>
+#include <grid/tria_accessor.h>
+#include <grid/dof_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary.h>
+#include <grid/dof_constraints.h>
+#include <basic/data_io.h>
+#include <basic/function.h>
+#include <basic/parameter_handler.h>
+#include <fe/fe_lib.h>
+#include <fe/quadrature_lib.h>
+#include <numerics/base.h>
+#include <numerics/assembler.h>
+#include <lac/dsmatrix.h>
+
+
+#include <map.h>
+#include <fstream.h>
+#include <cmath>
+extern "C" {
+# include <stdlib.h>
+}
+
+
+
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim> {
+ public:
+ /**
+ * Return the value of the function
+ * at the given point.
+ */
+ virtual double operator () (const Point<dim> &p) const;
+};
+
+
+
+
+
+
+template <int dim>
+class PoissonEquation : public Equation<dim> {
+ public:
+ PoissonEquation (const Function<dim> &rhs) :
+ Equation<dim>(1),
+ right_hand_side (rhs) {};
+
+ virtual void assemble (dFMatrix &cell_matrix,
+ dVector &rhs,
+ const FEValues<dim> &fe_values,
+ const Triangulation<dim>::cell_iterator &cell) const;
+ virtual void assemble (dFMatrix &cell_matrix,
+ const FEValues<dim> &fe_values,
+ const Triangulation<dim>::cell_iterator &cell) const;
+ virtual void assemble (dVector &rhs,
+ const FEValues<dim> &fe_values,
+ const Triangulation<dim>::cell_iterator &cell) const;
+ protected:
+ const Function<dim> &right_hand_side;
+};
+
+
+
+
+
+
+template <int dim>
+class PoissonProblem : public ProblemBase<dim>,
+ 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<dim> *tria;
+ DoFHandler<dim> *dof;
+};
+
+
+
+
+
+/*---------------------------- poisson.h ---------------------------*/
+/* end of #ifndef __poisson_H */
+#endif
+/*---------------------------- poisson.h ---------------------------*/
--- /dev/null
+/* $Id$ */
+
+
+#include "poisson.h"
+
+
+
+
+template <int dim>
+class CurvedLine :
+ public StraightBoundary<dim> {
+ public:
+ virtual Point<dim> in_between (const PointArray &neighbors) const {
+ Point<dim> middle = StraightBoundary<dim>::in_between(neighbors);
+ double x=middle(0),
+ y=middle(1);
+
+ if (y<x)
+ if (y<1-x)
+ middle(1) = 0.03*sin(6*3.141592*middle(0));
+ else
+ middle(0) = 1+0.03*sin(6*3.141592*middle(1));
+
+ else
+ if (y<1-x)
+ middle(0) = 0.03*sin(6*3.141592*middle(1));
+ else
+ middle(1) = 1+0.03*sin(6*3.141592*middle(0));
+
+ return middle;
+ };
+};
+
+
+
+
+template <int dim>
+PoissonProblem<dim>::PoissonProblem () :
+ tria(0), dof(0) {};
+
+
+
+
+template <int dim>
+void PoissonProblem<dim>::clear () {
+ if (tria != 0) delete tria;
+ if (dof != 0) delete dof;
+
+ ProblemBase<dim>::clear ();
+};
+
+
+
+
+template <int dim>
+void PoissonProblem<dim>::create_new (const unsigned int) {
+ clear ();
+
+ tria = new Triangulation<dim>();
+ dof = new DoFHandler<dim> (tria);
+ set_tria_and_dof (tria, dof);
+};
+
+
+
+
+template <int dim>
+void PoissonProblem<dim>::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 <int dim>
+bool PoissonProblem<dim>::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 <int dim>
+void PoissonProblem<dim>::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<dim>::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 <int dim>
+void PoissonProblem<dim>::make_ball_grid () {
+ // make ball grid around origin with
+ // unit radius
+ const Point<dim> origin;
+ static const HyperBallBoundary<dim> boundary(origin, 1.);
+
+ tria->create_hyper_ball (origin, 1.);
+ tria->set_boundary (&boundary);
+};
+
+
+
+
+template <int dim>
+void PoissonProblem<dim>::make_curved_line_grid () {
+ // set the boundary function
+ static const CurvedLine<dim> boundary;
+
+ tria->create_hypercube ();
+ tria->set_boundary (&boundary);
+};
+
+
+
+template <int dim>
+void PoissonProblem<dim>::make_random_grid () {
+ tria->refine_global (1);
+
+ Triangulation<dim>::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 <int dim>
+void PoissonProblem<dim>::run (ParameterHandler &prm) {
+ cout << "Test case = " << prm.get ("Test run") << endl
+ << endl;
+
+ cout << " Making grid..." << endl;
+ if (!make_grid (prm))
+ return;
+
+
+ FELinear<dim> fe;
+ static const RightHandSide<dim> rhs;
+ PoissonEquation<dim> equation (rhs);
+ QGauss4<dim> quadrature;
+
+ cout << "Distributing dofs... ";
+ dof->distribute_dofs (fe);
+ cout << dof->n_dofs() << " degrees of freedom." << endl;
+
+ cout << "Assembling matrices..." << endl;
+ FEValues<dim>::UpdateStruct update_flags;
+ update_flags.q_points = update_flags.gradients = true;
+ update_flags.jacobians = update_flags.JxW_values = true;
+
+ ProblemBase<dim>::DirichletBC dirichlet_bc;
+ ZeroFunction<dim> zero;
+ dirichlet_bc[0] = &zero;
+ assemble (equation, quadrature, fe, update_flags, dirichlet_bc);
+
+ cout << "Solving..." << endl;
+ solve ();
+
+ cout << "Printing..." << endl;
+ DataOut<dim> out;
+ ofstream gnuplot("gnuplot.out.5");
+ fill_data (out);
+ out.write_gnuplot (gnuplot);
+ gnuplot.close ();
+};
+
+
+
+
+
+template class PoissonProblem<1>;
+template class PoissonProblem<2>;
--- /dev/null
+# $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
--- /dev/null
+/* $Id$ */
+
+#include "poisson.h"
+
+
+
+template <int dim>
+double RightHandSide<dim>::operator () (const Point<dim> &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<fe_values.n_quadrature_points; ++point)
+ for (unsigned int i=0; i<fe_values.total_dofs; ++i)
+ {
+ for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+ cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
+ fe_values.shape_grad(j,point)) *
+ fe_values.JxW(point);
+ rhs(i) += fe_values.shape_value(i,point) *
+ right_hand_side(fe_values.quadrature_point(point)) *
+ fe_values.JxW(point);
+ };
+};
+
+
+
+void PoissonEquation<2>::assemble (dFMatrix &cell_matrix,
+ dVector &rhs,
+ const FEValues<2> &fe_values,
+ const Triangulation<2>::cell_iterator &) const {
+ for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+ for (unsigned int i=0; i<fe_values.total_dofs; ++i)
+ {
+ for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+ cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
+ fe_values.shape_grad(j,point)) *
+ fe_values.JxW(point);
+ rhs(i) += fe_values.shape_value(i,point) *
+ right_hand_side(fe_values.quadrature_point(point)) *
+ fe_values.JxW(point);
+ };
+};
+
+
+
+template <int dim>
+void PoissonEquation<dim>::assemble (dFMatrix &,
+ const FEValues<dim> &,
+ const Triangulation<dim>::cell_iterator &) const {
+ Assert (false, ExcPureVirtualFunctionCalled());
+};
+
+
+
+template <int dim>
+void PoissonEquation<dim>::assemble (dVector &,
+ const FEValues<dim> &,
+ const Triangulation<dim>::cell_iterator &) const {
+ Assert (false, ExcPureVirtualFunctionCalled());
+};
+
+
+
+
+
+template class RightHandSide<1>;
+template class RightHandSide<2>;
+
+
+template class PoissonEquation<1>;
+template class PoissonEquation<2>;
/* $Id$ */
-#include <grid/tria.h>
-#include <grid/dof.h>
-#include <grid/tria_accessor.h>
-#include <grid/dof_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_boundary.h>
-#include <grid/dof_constraints.h>
-#include <basic/data_io.h>
-#include <basic/function.h>
-#include <fe/fe_lib.h>
-#include <fe/quadrature_lib.h>
-#include <numerics/base.h>
-#include <numerics/assembler.h>
-#include <lac/dsmatrix.h>
+#include "poisson.h"
-#include <map.h>
-#include <fstream.h>
-#include <cmath>
-extern "C" {
-# include <stdlib.h>
-}
-extern TriaActiveIterator<1,CellAccessor<1> > x;
-extern TriaActiveIterator<2,CellAccessor<2> > y;
-extern TriaRawIterator<1,DoFLineAccessor<1,LineAccessor<1> > > z;
-
-
-
-template <int dim>
-class PoissonEquation : public Equation<dim> {
- public:
- PoissonEquation () :
- Equation<dim>(1) {};
-
- virtual void assemble (dFMatrix &cell_matrix,
- dVector &rhs,
- const FEValues<dim> &fe_values,
- const Triangulation<dim>::cell_iterator &cell) const;
- double right_hand_side (const Point<dim> &) const;
-};
-
-
-
-template <int dim>
-inline
-double PoissonEquation<dim>::right_hand_side (const Point<dim> &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<fe_values.n_quadrature_points; ++point)
- for (unsigned int i=0; i<fe_values.total_dofs; ++i)
- {
- for (unsigned int j=0; j<fe_values.total_dofs; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point) +
- fe_values.shape_value(i,point) *
- fe_values.shape_value(j,point)) *
- fe_values.JxW(point);
- rhs(i) += fe_values.shape_value(i,point) *
- right_hand_side(fe_values.quadrature_point(point)) *
- fe_values.JxW(point);
- };
-};
-
-
-
-void PoissonEquation<2>::assemble (dFMatrix &cell_matrix,
- dVector &rhs,
- const FEValues<2> &fe_values,
- const Triangulation<2>::cell_iterator &) const {
- for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
- for (unsigned int i=0; i<fe_values.total_dofs; ++i)
- {
- for (unsigned int j=0; j<fe_values.total_dofs; ++j)
- cell_matrix(i,j) += (fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point)/* +
- fe_values.shape_value(i,point) *
- fe_values.shape_value(j,point)*/) *
- fe_values.JxW(point);
- rhs(i) += fe_values.shape_value(i,point) *
- right_hand_side(fe_values.quadrature_point(point)) *
- fe_values.JxW(point);
- };
-};
-
-
-
-
-int main () {
- Triangulation<2> 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<dim>::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;
};
--- /dev/null
+/*---------------------------- poisson.h ---------------------------*/
+/* $Id$ */
+#ifndef __poisson_H
+#define __poisson_H
+/*---------------------------- poisson.h ---------------------------*/
+
+
+#include <grid/tria.h>
+#include <grid/dof.h>
+#include <grid/tria_accessor.h>
+#include <grid/dof_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary.h>
+#include <grid/dof_constraints.h>
+#include <basic/data_io.h>
+#include <basic/function.h>
+#include <basic/parameter_handler.h>
+#include <fe/fe_lib.h>
+#include <fe/quadrature_lib.h>
+#include <numerics/base.h>
+#include <numerics/assembler.h>
+#include <lac/dsmatrix.h>
+
+
+#include <map.h>
+#include <fstream.h>
+#include <cmath>
+extern "C" {
+# include <stdlib.h>
+}
+
+
+
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim> {
+ public:
+ /**
+ * Return the value of the function
+ * at the given point.
+ */
+ virtual double operator () (const Point<dim> &p) const;
+};
+
+
+
+
+
+
+template <int dim>
+class PoissonEquation : public Equation<dim> {
+ public:
+ PoissonEquation (const Function<dim> &rhs) :
+ Equation<dim>(1),
+ right_hand_side (rhs) {};
+
+ virtual void assemble (dFMatrix &cell_matrix,
+ dVector &rhs,
+ const FEValues<dim> &fe_values,
+ const Triangulation<dim>::cell_iterator &cell) const;
+ virtual void assemble (dFMatrix &cell_matrix,
+ const FEValues<dim> &fe_values,
+ const Triangulation<dim>::cell_iterator &cell) const;
+ virtual void assemble (dVector &rhs,
+ const FEValues<dim> &fe_values,
+ const Triangulation<dim>::cell_iterator &cell) const;
+ protected:
+ const Function<dim> &right_hand_side;
+};
+
+
+
+
+
+
+template <int dim>
+class PoissonProblem : public ProblemBase<dim>,
+ 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<dim> *tria;
+ DoFHandler<dim> *dof;
+};
+
+
+
+
+
+/*---------------------------- poisson.h ---------------------------*/
+/* end of #ifndef __poisson_H */
+#endif
+/*---------------------------- poisson.h ---------------------------*/
--- /dev/null
+/* $Id$ */
+
+
+#include "poisson.h"
+
+
+
+
+template <int dim>
+class CurvedLine :
+ public StraightBoundary<dim> {
+ public:
+ virtual Point<dim> in_between (const PointArray &neighbors) const {
+ Point<dim> middle = StraightBoundary<dim>::in_between(neighbors);
+ double x=middle(0),
+ y=middle(1);
+
+ if (y<x)
+ if (y<1-x)
+ middle(1) = 0.03*sin(6*3.141592*middle(0));
+ else
+ middle(0) = 1+0.03*sin(6*3.141592*middle(1));
+
+ else
+ if (y<1-x)
+ middle(0) = 0.03*sin(6*3.141592*middle(1));
+ else
+ middle(1) = 1+0.03*sin(6*3.141592*middle(0));
+
+ return middle;
+ };
+};
+
+
+
+
+template <int dim>
+PoissonProblem<dim>::PoissonProblem () :
+ tria(0), dof(0) {};
+
+
+
+
+template <int dim>
+void PoissonProblem<dim>::clear () {
+ if (tria != 0) delete tria;
+ if (dof != 0) delete dof;
+
+ ProblemBase<dim>::clear ();
+};
+
+
+
+
+template <int dim>
+void PoissonProblem<dim>::create_new (const unsigned int) {
+ clear ();
+
+ tria = new Triangulation<dim>();
+ dof = new DoFHandler<dim> (tria);
+ set_tria_and_dof (tria, dof);
+};
+
+
+
+
+template <int dim>
+void PoissonProblem<dim>::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 <int dim>
+bool PoissonProblem<dim>::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 <int dim>
+void PoissonProblem<dim>::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<dim>::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 <int dim>
+void PoissonProblem<dim>::make_ball_grid () {
+ // make ball grid around origin with
+ // unit radius
+ const Point<dim> origin;
+ static const HyperBallBoundary<dim> boundary(origin, 1.);
+
+ tria->create_hyper_ball (origin, 1.);
+ tria->set_boundary (&boundary);
+};
+
+
+
+
+template <int dim>
+void PoissonProblem<dim>::make_curved_line_grid () {
+ // set the boundary function
+ static const CurvedLine<dim> boundary;
+
+ tria->create_hypercube ();
+ tria->set_boundary (&boundary);
+};
+
+
+
+template <int dim>
+void PoissonProblem<dim>::make_random_grid () {
+ tria->refine_global (1);
+
+ Triangulation<dim>::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 <int dim>
+void PoissonProblem<dim>::run (ParameterHandler &prm) {
+ cout << "Test case = " << prm.get ("Test run") << endl
+ << endl;
+
+ cout << " Making grid..." << endl;
+ if (!make_grid (prm))
+ return;
+
+
+ FELinear<dim> fe;
+ static const RightHandSide<dim> rhs;
+ PoissonEquation<dim> equation (rhs);
+ QGauss4<dim> quadrature;
+
+ cout << "Distributing dofs... ";
+ dof->distribute_dofs (fe);
+ cout << dof->n_dofs() << " degrees of freedom." << endl;
+
+ cout << "Assembling matrices..." << endl;
+ FEValues<dim>::UpdateStruct update_flags;
+ update_flags.q_points = update_flags.gradients = true;
+ update_flags.jacobians = update_flags.JxW_values = true;
+
+ ProblemBase<dim>::DirichletBC dirichlet_bc;
+ ZeroFunction<dim> zero;
+ dirichlet_bc[0] = &zero;
+ assemble (equation, quadrature, fe, update_flags, dirichlet_bc);
+
+ cout << "Solving..." << endl;
+ solve ();
+
+ cout << "Printing..." << endl;
+ DataOut<dim> out;
+ ofstream gnuplot("gnuplot.out.5");
+ fill_data (out);
+ out.write_gnuplot (gnuplot);
+ gnuplot.close ();
+};
+
+
+
+
+
+template class PoissonProblem<1>;
+template class PoissonProblem<2>;