# $Id$
CXX=c++
-INCPATH = -I../include -I../../lac/include -I../../base/include
-CXXFLAGS.g= -DDEBUG -g -Wall $(INCPATH)
-CXXFLAGS = -O3 $(INCPATH)
+INCLUDE = -I../include -I../../lac/include -I../../base/include
+CXXFLAGS.g= -DDEBUG -g -Wall -W -pedantic -Wconversion \
+ -Winline -Woverloaded-virtual -fno-rtti -fno-exceptions \
+ $(INCLUDE)
+CXXFLAGS =-O3 -Wuninitialized -finline-functions -ffast-math \
+ -funroll-loops -felide-constructors -fnonnull-objects \
+ -fno-rtti -fno-exceptions $(INCLUDE)
LIBPATH = -L../lib -L../../lac/lib
LIBFILES.g= ../lib/libbasic.g.a ../lib/libgrid.g.a \
LIBS.g = $(LIBPATH) -lbasic.g -lgrid.g -lfe.g -lnumerics.g -llac.g
LIBS = $(LIBS.g:.g=)
+
+ifeq ($(shell uname),Linux)
+CXX = /home/wolf/bin/gcc-2.8.1/bin/c++
+endif
+
+
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)
+CXXFLAGS := $(CXXFLAGS) $(INCLUDE2.8)
+CXXFLAGS.g:= $(CXXFLAGS.g) $(INCLUDE2.8)
endif
#of this file.
Makefile.dep: $(cc-files) $(h-files)
@echo ================= Making $@
- @$(CXX) $(INCPATH) -M $(cc-files) > Makefile.dep
+ @$(CXX) $(INCLUDE) -M $(cc-files) > Makefile.dep
include Makefile.dep
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.g= -DDEBUG -g -Wall -W -pedantic -Wconversion \
+ -Winline -Woverloaded-virtual -fno-rtti -fno-exceptions \
+ $(INCLUDE)
CXXFLAGS =-O3 -Wuninitialized -finline-functions -ffast-math \
- -funroll-loops -felide-constructors -fnonnull-objects $(INCLUDE)
+ -funroll-loops -felide-constructors -fnonnull-objects \
+ -fno-rtti -fno-exceptions $(INCLUDE)
+
+ifeq ($(shell uname),Linux)
+CXX = /home/wolf/bin/gcc-2.8.1/bin/c++
+endif
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)
+CXXFLAGS := $(CXXFLAGS) $(INCLUDE2.8)
+CXXFLAGS.g:= $(CXXFLAGS.g) $(INCLUDE2.8)
endif
LIBPATH = -L../../lib -L../../../lac/lib
-set Test run = { zoom in | ball | curved line | random }
-set Global refinement = {{ 2 | 5 | 6 | 0 }}
-set Right hand side = {{ zero | zero | trigpoly | constant }}
-set Boundary values = {{ sine | sine | zero | zero }}
-set Output file = gnuplot.{{ zoom_in | ball | curved_line | random }}
+set Test run = { zoom in | ball | curved line | random | jump | L-region}
+set Global refinement = {{ 2 | 5 | 6 | 0 | 2 | 5 }}
+set Right hand side = {{ zero | zero | trigpoly | constant | zero | zero }}
+set Boundary values = {{ sine | sine | zero | zero | jump | sine }}
+set Output file = gnuplot.{{ zoom_in | ball | curved_line | random | jump | L-region }}
+template <int dim>
+class BoundaryValuesJump : public Function<dim> {
+ public:
+ /**
+ * Return the value of the function
+ * at the given point.
+ */
+ virtual double operator () (const Point<dim> &p) const {
+ switch (dim)
+ {
+ case 1:
+ return 0;
+ default:
+ if (p(0) == p(1))
+ return 0.5;
+ else
+ return (p(0)>p(1) ? 0. : 1.);
+ };
+ };
+};
+
+
+
template <int dim>
class RHSTrigPoly : public Function<dim> {
+/**
+ Right hand side constructed such that the exact solution is
+ $x(1-x)$ in 1d, $x(1-x)*y(1-y)$ in 2d, etc.
+ */
+template <int dim>
+class RHSPoly : 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 Solution : 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>
+double RHSPoly<dim>::operator () (const Point<dim> &p) const {
+ double ret_val = 0;
+ for (unsigned int i=0; i<dim; ++i)
+ ret_val += 2*p(i)*(1.-p(i));
+ return ret_val;
+};
+
+
+template <int dim>
+double Solution<dim>::operator () (const Point<dim> &p) const {
+ double ret_val = 0;
+ for (unsigned int i=0; i<dim; ++i)
+ ret_val *= p(i)*(1.-p(i));
+ return ret_val;
+};
+
+
+
template <int dim>
void PoissonProblem<dim>::declare_parameters (ParameterHandler &prm) {
if (dim>=2)
prm.declare_entry ("Test run", "zoom in",
- "tensor\\|zoom in\\|ball\\|curved line\\|random");
+ "tensor\\|zoom in\\|ball\\|curved line\\|random\\|jump\\|L-region");
else
prm.declare_entry ("Test run", "zoom in", "tensor\\|zoom in\\|random");
prm.declare_entry ("Global refinement", "0",
ParameterHandler::RegularExpressions::Integer);
prm.declare_entry ("Right hand side", "zero",
- "zero\\|constant\\|trigpoly");
+ "zero\\|constant\\|trigpoly\\|poly");
prm.declare_entry ("Boundary values", "zero",
- "zero\\|sine");
+ "zero\\|sine\\|jump");
prm.declare_entry ("Output file", "gnuplot.1");
};
else
if (test=="tensor") test_case = 5;
else
- {
- cerr << "This test seems not to be implemented!" << endl;
- return false;
- };
+ if (test=="jump") test_case = 6;
+ else
+ if (test=="L-region") test_case = 7;
+ else
+ {
+ cerr << "This test seems not to be implemented!" << endl;
+ return false;
+ };
switch (test_case)
{
case 5:
tria->create_hypercube ();
break;
+ case 6:
+ tria->create_hypercube ();
+ tria->refine_global (1);
+ for (unsigned int i=0; i<5; ++i)
+ {
+ tria->begin_active(tria->n_levels()-1)->set_refine_flag();
+ (--(tria->last_active()))->set_refine_flag();
+ tria->execute_refinement ();
+ };
+ break;
+ case 7:
+ tria->create_hyper_L ();
+ break;
default:
return false;
};
if (rhs_name == "trigpoly")
rhs = new RHSTrigPoly<dim>();
else
- return false;
+ if (rhs_name == "poly")
+ rhs = new RHSPoly<dim> ();
+ else
+ return false;
if (rhs != 0)
return true;
if (bv_name == "sine")
boundary_values = new BoundaryValuesSine<dim> ();
else
- return false;
+ if (bv_name == "jump")
+ boundary_values = new BoundaryValuesJump<dim> ();
+ else
+ return false;
if (boundary_values != 0)
return true;
cout << "Test case = " << prm.get ("Test run")
<< endl;
- cout << " Making grid..." << endl;
+ cout << " Making grid... ";
if (!make_grid (prm))
return;
+ cout << tria->n_active_cells() << " active cells." << endl;
if (!set_right_hand_side (prm))
return;
FELinear<dim> fe;
PoissonEquation<dim> equation (*rhs);
- QGauss4<dim> quadrature;
+ QGauss3<dim> quadrature;
cout << " Distributing dofs... ";
dof->distribute_dofs (fe);
solve ();
cout << " Writing to file <" << prm.get("Output file") << ">..."
- << endl
<< endl;
DataOut<dim> out;
- ofstream gnuplot(prm.get("Output file"));
+ String o_filename = prm.get ("Output file");
+ ofstream gnuplot(o_filename);
fill_data (out);
out.write_gnuplot (gnuplot);
gnuplot.close ();
+
+ cout << " Calculation L2 error... ";
+ Solution<dim> sol;
+ vector<double> difference;
+ QGauss3<dim> q;
+ integrate_difference (sol, difference, q, fe, L2_norm);
+ double error = 0;
+ for (unsigned int i=0; i<difference.size(); ++i)
+ error += difference[i]*difference[i];
+ cout << sqrt(error) << endl;
+
+ cout << endl;
};
# $Id$
CXX=c++
-INCPATH = -I../include -I../../lac/include -I../../base/include
-CXXFLAGS.g= -DDEBUG -g -Wall $(INCPATH)
-CXXFLAGS = -O3 $(INCPATH)
+INCLUDE = -I../include -I../../lac/include -I../../base/include
+CXXFLAGS.g= -DDEBUG -g -Wall -W -pedantic -Wconversion \
+ -Winline -Woverloaded-virtual -fno-rtti -fno-exceptions \
+ $(INCLUDE)
+CXXFLAGS =-O3 -Wuninitialized -finline-functions -ffast-math \
+ -funroll-loops -felide-constructors -fnonnull-objects \
+ -fno-rtti -fno-exceptions $(INCLUDE)
LIBPATH = -L../lib -L../../lac/lib
LIBFILES.g= ../lib/libbasic.g.a ../lib/libgrid.g.a \
LIBS.g = $(LIBPATH) -lbasic.g -lgrid.g -lfe.g -lnumerics.g -llac.g
LIBS = $(LIBS.g:.g=)
+
+ifeq ($(shell uname),Linux)
+CXX = /home/wolf/bin/gcc-2.8.1/bin/c++
+endif
+
+
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)
+CXXFLAGS := $(CXXFLAGS) $(INCLUDE2.8)
+CXXFLAGS.g:= $(CXXFLAGS.g) $(INCLUDE2.8)
endif
#of this file.
Makefile.dep: $(cc-files) $(h-files)
@echo ================= Making $@
- @$(CXX) $(INCPATH) -M $(cc-files) > Makefile.dep
+ @$(CXX) $(INCLUDE) -M $(cc-files) > Makefile.dep
include Makefile.dep
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.g= -DDEBUG -g -Wall -W -pedantic -Wconversion \
+ -Winline -Woverloaded-virtual -fno-rtti -fno-exceptions \
+ $(INCLUDE)
CXXFLAGS =-O3 -Wuninitialized -finline-functions -ffast-math \
- -funroll-loops -felide-constructors -fnonnull-objects $(INCLUDE)
+ -funroll-loops -felide-constructors -fnonnull-objects \
+ -fno-rtti -fno-exceptions $(INCLUDE)
+
+ifeq ($(shell uname),Linux)
+CXX = /home/wolf/bin/gcc-2.8.1/bin/c++
+endif
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)
+CXXFLAGS := $(CXXFLAGS) $(INCLUDE2.8)
+CXXFLAGS.g:= $(CXXFLAGS.g) $(INCLUDE2.8)
endif
LIBPATH = -L../../lib -L../../../lac/lib
set output "fig.random.eps"
splot "gnuplot.random"
+!echo " Making <fig.jump.eps>"
+set output "fig.jump.eps"
+splot "gnuplot.jump"
+!echo " Making <fig.L-region.eps>"
+set view 52,115
+set output "fig.L-region.eps"
+splot "gnuplot.L-region"
-set Test run = { zoom in | ball | curved line | random }
-set Global refinement = {{ 2 | 5 | 6 | 0 }}
-set Right hand side = {{ zero | zero | trigpoly | constant }}
-set Boundary values = {{ sine | sine | zero | zero }}
-set Output file = gnuplot.{{ zoom_in | ball | curved_line | random }}
+set Test run = { zoom in | ball | curved line | random | jump | L-region}
+set Global refinement = {{ 2 | 5 | 6 | 0 | 2 | 5 }}
+set Right hand side = {{ zero | zero | trigpoly | constant | zero | zero }}
+set Boundary values = {{ sine | sine | zero | zero | jump | sine }}
+set Output file = gnuplot.{{ zoom_in | ball | curved_line | random | jump | L-region }}
+template <int dim>
+class BoundaryValuesJump : public Function<dim> {
+ public:
+ /**
+ * Return the value of the function
+ * at the given point.
+ */
+ virtual double operator () (const Point<dim> &p) const {
+ switch (dim)
+ {
+ case 1:
+ return 0;
+ default:
+ if (p(0) == p(1))
+ return 0.5;
+ else
+ return (p(0)>p(1) ? 0. : 1.);
+ };
+ };
+};
+
+
+
template <int dim>
class RHSTrigPoly : public Function<dim> {
+/**
+ Right hand side constructed such that the exact solution is
+ $x(1-x)$ in 1d, $x(1-x)*y(1-y)$ in 2d, etc.
+ */
+template <int dim>
+class RHSPoly : 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 Solution : 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>
+double RHSPoly<dim>::operator () (const Point<dim> &p) const {
+ double ret_val = 0;
+ for (unsigned int i=0; i<dim; ++i)
+ ret_val += 2*p(i)*(1.-p(i));
+ return ret_val;
+};
+
+
+template <int dim>
+double Solution<dim>::operator () (const Point<dim> &p) const {
+ double ret_val = 0;
+ for (unsigned int i=0; i<dim; ++i)
+ ret_val *= p(i)*(1.-p(i));
+ return ret_val;
+};
+
+
+
template <int dim>
void PoissonProblem<dim>::declare_parameters (ParameterHandler &prm) {
if (dim>=2)
prm.declare_entry ("Test run", "zoom in",
- "tensor\\|zoom in\\|ball\\|curved line\\|random");
+ "tensor\\|zoom in\\|ball\\|curved line\\|random\\|jump\\|L-region");
else
prm.declare_entry ("Test run", "zoom in", "tensor\\|zoom in\\|random");
prm.declare_entry ("Global refinement", "0",
ParameterHandler::RegularExpressions::Integer);
prm.declare_entry ("Right hand side", "zero",
- "zero\\|constant\\|trigpoly");
+ "zero\\|constant\\|trigpoly\\|poly");
prm.declare_entry ("Boundary values", "zero",
- "zero\\|sine");
+ "zero\\|sine\\|jump");
prm.declare_entry ("Output file", "gnuplot.1");
};
else
if (test=="tensor") test_case = 5;
else
- {
- cerr << "This test seems not to be implemented!" << endl;
- return false;
- };
+ if (test=="jump") test_case = 6;
+ else
+ if (test=="L-region") test_case = 7;
+ else
+ {
+ cerr << "This test seems not to be implemented!" << endl;
+ return false;
+ };
switch (test_case)
{
case 5:
tria->create_hypercube ();
break;
+ case 6:
+ tria->create_hypercube ();
+ tria->refine_global (1);
+ for (unsigned int i=0; i<5; ++i)
+ {
+ tria->begin_active(tria->n_levels()-1)->set_refine_flag();
+ (--(tria->last_active()))->set_refine_flag();
+ tria->execute_refinement ();
+ };
+ break;
+ case 7:
+ tria->create_hyper_L ();
+ break;
default:
return false;
};
if (rhs_name == "trigpoly")
rhs = new RHSTrigPoly<dim>();
else
- return false;
+ if (rhs_name == "poly")
+ rhs = new RHSPoly<dim> ();
+ else
+ return false;
if (rhs != 0)
return true;
if (bv_name == "sine")
boundary_values = new BoundaryValuesSine<dim> ();
else
- return false;
+ if (bv_name == "jump")
+ boundary_values = new BoundaryValuesJump<dim> ();
+ else
+ return false;
if (boundary_values != 0)
return true;
cout << "Test case = " << prm.get ("Test run")
<< endl;
- cout << " Making grid..." << endl;
+ cout << " Making grid... ";
if (!make_grid (prm))
return;
+ cout << tria->n_active_cells() << " active cells." << endl;
if (!set_right_hand_side (prm))
return;
FELinear<dim> fe;
PoissonEquation<dim> equation (*rhs);
- QGauss4<dim> quadrature;
+ QGauss3<dim> quadrature;
cout << " Distributing dofs... ";
dof->distribute_dofs (fe);
solve ();
cout << " Writing to file <" << prm.get("Output file") << ">..."
- << endl
<< endl;
DataOut<dim> out;
- ofstream gnuplot(prm.get("Output file"));
+ String o_filename = prm.get ("Output file");
+ ofstream gnuplot(o_filename);
fill_data (out);
out.write_gnuplot (gnuplot);
gnuplot.close ();
+
+ cout << " Calculation L2 error... ";
+ Solution<dim> sol;
+ vector<double> difference;
+ QGauss3<dim> q;
+ integrate_difference (sol, difference, q, fe, L2_norm);
+ double error = 0;
+ for (unsigned int i=0; i<difference.size(); ++i)
+ error += difference[i]*difference[i];
+ cout << sqrt(error) << endl;
+
+ cout << endl;
};