]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add an example for nonlinear problems.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 25 Feb 1999 10:26:39 +0000 (10:26 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 25 Feb 1999 10:26:39 +0000 (10:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@905 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/Attic/examples/nonlinear/fixed-point-iteration/Makefile [new file with mode: 0644]
deal.II/deal.II/Attic/examples/nonlinear/fixed-point-iteration/nonlinear.cc [new file with mode: 0644]
tests/big-tests/nonlinear/fixed-point-iteration/Makefile [new file with mode: 0644]
tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc [new file with mode: 0644]

diff --git a/deal.II/deal.II/Attic/examples/nonlinear/fixed-point-iteration/Makefile b/deal.II/deal.II/Attic/examples/nonlinear/fixed-point-iteration/Makefile
new file mode 100644 (file)
index 0000000..949f8f5
--- /dev/null
@@ -0,0 +1,122 @@
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1998
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly:
+target   = nonlinear
+
+# All dependencies between files should be updated by the included
+# file Makefile.dep if necessary. Object files are compiled into
+# the archives ./Obj.a and ./Obj.g.a. By default, the debug version
+# is used to link. It you don't like that, change the following
+# variable to "off"
+debug-mode = on
+
+# If you want your program to be linked with extra object or library
+# files, specify them here:
+user-libs =
+
+# To run the program, use "make run"; to give parameters to the program,
+# give the parameters to the following variable:
+run-parameters  = $(target).prm
+
+# To execute additional action apart from running the program, fill
+# in this list:
+additional-run-action = gnuplot make_ps
+
+# To specify which files are to be deleted by "make clean" (apart from
+# the usual ones: object files, executables, backups, etc), fill in the
+# following list
+delete-files = gnuplot* *.eps
+
+
+
+
+###############################################################################
+# Internals
+
+#deal include base path
+D = ../../../..
+
+include ../../../Make.global_options
+
+
+
+# get lists of files we need
+cc-files    = $(wildcard *.cc)
+o-files     = $(cc-files:.cc=.o)
+go-files    = $(cc-files:.cc=.go)
+h-files     = $(wildcard *.h)
+lib-h-files = $(wildcard ../../include/*/*.h)
+
+# list of libraries needed to link with
+libs     = ./Obj.a   -ldeal_II_2d   -llac -lbase
+libs.g   = ./Obj.g.a -ldeal_II_2d.g -llac.g -lbase.g
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+libraries = $(libs.g)
+flags     = $(CXXFLAGS.g)
+endif
+
+ifeq ($(debug-mode),off)
+libraries = $(libs)
+flags     = $(CXXFLAGS)
+endif
+
+
+
+# make rule for the target
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) $(flags) -o $@ $^ $(user-libs)
+
+# rule how to run the program
+run: $(target)
+       $(target) $(run-parameters)
+       $(additional-run-action)
+
+
+# rule to make object files
+%.go : %.cc
+       @echo ============================ Compiling with debugging information:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+       @echo ============================ Compiling with optimization:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+# rules which files the libraries depend upon
+Obj.a: ./Obj.a($(o-files))
+Obj.g.a: ./Obj.g.a($(go-files))
+
+
+clean:
+       -rm -f *.o *.go *~ Makefile.dep Obj.a Obj.g.a $(target) $(delete-files)
+
+
+
+.PHONY: clean
+
+
+#Rule to generate the dependency file. This file is
+#automagically remade whenever needed, i.e. whenever
+#one of the cc-/h-files changed. Make detects whether
+#to remake this file upon inclusion at the bottom
+#of this file.
+#
+#use perl to generate rules for the .go files as well
+#as to make rules not for tria.o and the like, but
+#rather for libnumerics.a(tria.o)
+Makefile.dep: $(cc-files) $(h-files) $(lib-h-files)
+       @echo ============================ Remaking Makefile
+       @perl ../../../Make_dep.pl ./Obj $(INCLUDE) $(cc-files) \
+               > Makefile.dep
+
+
+include Makefile.dep
+
diff --git a/deal.II/deal.II/Attic/examples/nonlinear/fixed-point-iteration/nonlinear.cc b/deal.II/deal.II/Attic/examples/nonlinear/fixed-point-iteration/nonlinear.cc
new file mode 100644 (file)
index 0000000..62bec30
--- /dev/null
@@ -0,0 +1,232 @@
+/* $Id$ */
+/* Copyright W. Bangerth, University of Heidelberg, 1998 */
+
+
+
+
+#include <lac/vector.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 <basic/data_io.h>
+#include <base/function.h>
+#include <fe/fe_lib.lagrange.h>
+#include <base/quadrature_lib.h>
+#include <numerics/base.h>
+#include <numerics/assembler.h>
+
+#include <map>
+#include <fstream>
+#include <cmath>
+#include <cstdlib>
+
+
+
+
+
+template <int dim>
+class RightHandSide :  public Function<dim> 
+{
+  public:
+    double operator () (const Point<dim> &p) const 
+      {
+       double x = 6;
+       for (unsigned int d=0; d<dim; ++d)
+         if (p(d) < 0.5)
+           x *= -1;
+       
+       return x;
+      };
+};
+
+
+
+template <int dim>
+class PoissonEquation :  public Equation<dim> {
+  public:
+    PoissonEquation (const Function<dim>  &rhs,
+                    const Vector<double> &last_solution) :
+                   Equation<dim>(1),
+                   right_hand_side (rhs),
+                   last_solution(last_solution) {};
+
+    virtual void assemble (FullMatrix<double>  &cell_matrix,
+                          Vector<double>      &rhs,
+                          const FEValues<dim> &fe_values,
+                          const DoFHandler<dim>::cell_iterator &cell) const;
+    virtual void assemble (FullMatrix<double>  &cell_matrix,
+                          const FEValues<dim> &fe_values,
+                          const DoFHandler<dim>::cell_iterator &cell) const;
+    virtual void assemble (Vector<double>      &rhs,
+                          const FEValues<dim> &fe_values,
+                          const DoFHandler<dim>::cell_iterator &cell) const;
+  protected:
+    const Function<dim>  &right_hand_side;
+    const Vector<double> &last_solution;
+};
+
+
+
+
+
+
+template <int dim>
+class NonlinearProblem : public ProblemBase<dim> {
+  public:
+    NonlinearProblem ();
+    void run ();
+
+  protected:
+    Triangulation<dim> *tria;
+    DoFHandler<dim>    *dof;
+
+    Vector<double>      last_solution;
+};
+
+
+
+
+template <int dim>
+void PoissonEquation<dim>::assemble (FullMatrix<double>  &cell_matrix,
+                                    Vector<double>      &rhs,
+                                    const FEValues<dim> &fe_values,
+                                    const DoFHandler<dim>::cell_iterator &) const {
+  const vector<vector<Tensor<1,dim> > >&gradients = fe_values.get_shape_grads ();
+  const FullMatrix<double> &values    = fe_values.get_shape_values ();
+  vector<double>            rhs_values (fe_values.n_quadrature_points);
+  const vector<double> &weights   = fe_values.get_JxW_values ();
+  
+  vector<Tensor<1,dim> >   last_solution_grads(fe_values.n_quadrature_points);
+  fe_values.get_function_grads (last_solution, last_solution_grads);
+
+  
+  right_hand_side.value_list (fe_values.get_quadrature_points(), rhs_values);
+   
+  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) += (gradients[i][point] *
+                              gradients[j][point]) *
+                             weights[point] /
+                             sqrt(1+last_solution_grads[i]*last_solution_grads[i]);
+       rhs(i) += values(i,point) *
+                 rhs_values[point] *
+                 weights[point];
+      };
+};
+
+
+
+
+template <int dim>
+void PoissonEquation<dim>::assemble (FullMatrix<double>  &,
+                                    const FEValues<dim> &,
+                                    const DoFHandler<dim>::cell_iterator &) const {
+  Assert (false, ExcPureVirtualFunctionCalled());
+};
+
+
+
+template <int dim>
+void PoissonEquation<dim>::assemble (Vector<double>      &,
+                                    const FEValues<dim> &,
+                                    const DoFHandler<dim>::cell_iterator &) const {
+  Assert (false, ExcPureVirtualFunctionCalled());
+};
+
+
+
+
+template <int dim>
+NonlinearProblem<dim>::NonlinearProblem () :
+               tria(0), dof(0) {};
+
+
+
+template <int dim>
+void NonlinearProblem<dim>::run () {
+
+                                  // first reset everything to a virgin state
+  clear ();
+  
+  tria = new Triangulation<dim>();
+  dof = new DoFHandler<dim> (tria);
+  set_tria_and_dof (tria, dof);
+
+
+  RightHandSide<dim>    rhs;
+  ZeroFunction<dim>     boundary_values;
+  StraightBoundary<dim> boundary;
+  
+  FELinear<dim>                   fe;
+  PoissonEquation<dim>            equation (rhs, last_solution);
+  QGauss2<dim>                    quadrature;
+  
+  ProblemBase<dim>::FunctionMap dirichlet_bc;
+  dirichlet_bc[0] = &boundary_values;
+
+
+  cout << "    Making grid... ";
+  tria->create_hypercube ();
+  tria->refine_global (6);
+  cout << tria->n_active_cells() << " active cells." << endl;
+  
+  cout << "    Distributing dofs... "; 
+  dof->distribute_dofs (fe);
+  cout << dof->n_dofs() << " degrees of freedom." << endl;
+
+                                  // set the starting values for the iteration
+                                  // to a constant value of 1
+  last_solution.reinit (dof->n_dofs());
+  for (unsigned int i=0; i<dof->n_dofs(); ++i)
+    last_solution(i) = 1;
+  
+
+                                  // here comes the fixed point iteration
+  for (unsigned int nonlinear_step=0; nonlinear_step<40; ++nonlinear_step)
+    {
+      cout << "    Nonlinear step " << nonlinear_step << endl;
+      cout << "        Assembling matrices..." << endl;
+      assemble (equation, quadrature, fe,
+               UpdateFlags(update_gradients | update_JxW_values |
+                           update_q_points),
+               dirichlet_bc, boundary);
+
+      cout << "        Solving..." << endl;
+      solve ();
+
+      if (nonlinear_step % 4 == 0)
+       {
+         string filename = "nonlinear.";
+         filename += ('0' + (nonlinear_step/4));
+         filename += ".gnuplot";
+         cout << "        Writing to file <" << filename << ">..." << endl;
+         
+         DataOut<dim> out;
+         ofstream gnuplot(filename.c_str());
+         fill_data (out);
+         out.write_gnuplot (gnuplot);
+         gnuplot.close ();
+       };
+
+      last_solution = solution;
+    };
+  
+  delete dof;
+  delete tria;
+  
+  cout << endl;
+};
+
+
+
+
+int main ()
+{
+  NonlinearProblem<2> problem;
+  problem.run ();
+};
diff --git a/tests/big-tests/nonlinear/fixed-point-iteration/Makefile b/tests/big-tests/nonlinear/fixed-point-iteration/Makefile
new file mode 100644 (file)
index 0000000..949f8f5
--- /dev/null
@@ -0,0 +1,122 @@
+# $Id$
+# Copyright W. Bangerth, University of Heidelberg, 1998
+
+# Template for makefiles for the examples subdirectory. In principle,
+# everything should be done automatically if you set the target file
+# here correctly:
+target   = nonlinear
+
+# All dependencies between files should be updated by the included
+# file Makefile.dep if necessary. Object files are compiled into
+# the archives ./Obj.a and ./Obj.g.a. By default, the debug version
+# is used to link. It you don't like that, change the following
+# variable to "off"
+debug-mode = on
+
+# If you want your program to be linked with extra object or library
+# files, specify them here:
+user-libs =
+
+# To run the program, use "make run"; to give parameters to the program,
+# give the parameters to the following variable:
+run-parameters  = $(target).prm
+
+# To execute additional action apart from running the program, fill
+# in this list:
+additional-run-action = gnuplot make_ps
+
+# To specify which files are to be deleted by "make clean" (apart from
+# the usual ones: object files, executables, backups, etc), fill in the
+# following list
+delete-files = gnuplot* *.eps
+
+
+
+
+###############################################################################
+# Internals
+
+#deal include base path
+D = ../../../..
+
+include ../../../Make.global_options
+
+
+
+# get lists of files we need
+cc-files    = $(wildcard *.cc)
+o-files     = $(cc-files:.cc=.o)
+go-files    = $(cc-files:.cc=.go)
+h-files     = $(wildcard *.h)
+lib-h-files = $(wildcard ../../include/*/*.h)
+
+# list of libraries needed to link with
+libs     = ./Obj.a   -ldeal_II_2d   -llac -lbase
+libs.g   = ./Obj.g.a -ldeal_II_2d.g -llac.g -lbase.g
+
+
+# check whether we use debug mode or not
+ifeq ($(debug-mode),on)
+libraries = $(libs.g)
+flags     = $(CXXFLAGS.g)
+endif
+
+ifeq ($(debug-mode),off)
+libraries = $(libs)
+flags     = $(CXXFLAGS)
+endif
+
+
+
+# make rule for the target
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) $(flags) -o $@ $^ $(user-libs)
+
+# rule how to run the program
+run: $(target)
+       $(target) $(run-parameters)
+       $(additional-run-action)
+
+
+# rule to make object files
+%.go : %.cc
+       @echo ============================ Compiling with debugging information:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc
+       @echo ============================ Compiling with optimization:   $<
+       @echo $(CXX) ... -c $< -o $@
+       @$(CXX) $(CXXFLAGS) -c $< -o $@
+
+
+# rules which files the libraries depend upon
+Obj.a: ./Obj.a($(o-files))
+Obj.g.a: ./Obj.g.a($(go-files))
+
+
+clean:
+       -rm -f *.o *.go *~ Makefile.dep Obj.a Obj.g.a $(target) $(delete-files)
+
+
+
+.PHONY: clean
+
+
+#Rule to generate the dependency file. This file is
+#automagically remade whenever needed, i.e. whenever
+#one of the cc-/h-files changed. Make detects whether
+#to remake this file upon inclusion at the bottom
+#of this file.
+#
+#use perl to generate rules for the .go files as well
+#as to make rules not for tria.o and the like, but
+#rather for libnumerics.a(tria.o)
+Makefile.dep: $(cc-files) $(h-files) $(lib-h-files)
+       @echo ============================ Remaking Makefile
+       @perl ../../../Make_dep.pl ./Obj $(INCLUDE) $(cc-files) \
+               > Makefile.dep
+
+
+include Makefile.dep
+
diff --git a/tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc b/tests/big-tests/nonlinear/fixed-point-iteration/nonlinear.cc
new file mode 100644 (file)
index 0000000..62bec30
--- /dev/null
@@ -0,0 +1,232 @@
+/* $Id$ */
+/* Copyright W. Bangerth, University of Heidelberg, 1998 */
+
+
+
+
+#include <lac/vector.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 <basic/data_io.h>
+#include <base/function.h>
+#include <fe/fe_lib.lagrange.h>
+#include <base/quadrature_lib.h>
+#include <numerics/base.h>
+#include <numerics/assembler.h>
+
+#include <map>
+#include <fstream>
+#include <cmath>
+#include <cstdlib>
+
+
+
+
+
+template <int dim>
+class RightHandSide :  public Function<dim> 
+{
+  public:
+    double operator () (const Point<dim> &p) const 
+      {
+       double x = 6;
+       for (unsigned int d=0; d<dim; ++d)
+         if (p(d) < 0.5)
+           x *= -1;
+       
+       return x;
+      };
+};
+
+
+
+template <int dim>
+class PoissonEquation :  public Equation<dim> {
+  public:
+    PoissonEquation (const Function<dim>  &rhs,
+                    const Vector<double> &last_solution) :
+                   Equation<dim>(1),
+                   right_hand_side (rhs),
+                   last_solution(last_solution) {};
+
+    virtual void assemble (FullMatrix<double>  &cell_matrix,
+                          Vector<double>      &rhs,
+                          const FEValues<dim> &fe_values,
+                          const DoFHandler<dim>::cell_iterator &cell) const;
+    virtual void assemble (FullMatrix<double>  &cell_matrix,
+                          const FEValues<dim> &fe_values,
+                          const DoFHandler<dim>::cell_iterator &cell) const;
+    virtual void assemble (Vector<double>      &rhs,
+                          const FEValues<dim> &fe_values,
+                          const DoFHandler<dim>::cell_iterator &cell) const;
+  protected:
+    const Function<dim>  &right_hand_side;
+    const Vector<double> &last_solution;
+};
+
+
+
+
+
+
+template <int dim>
+class NonlinearProblem : public ProblemBase<dim> {
+  public:
+    NonlinearProblem ();
+    void run ();
+
+  protected:
+    Triangulation<dim> *tria;
+    DoFHandler<dim>    *dof;
+
+    Vector<double>      last_solution;
+};
+
+
+
+
+template <int dim>
+void PoissonEquation<dim>::assemble (FullMatrix<double>  &cell_matrix,
+                                    Vector<double>      &rhs,
+                                    const FEValues<dim> &fe_values,
+                                    const DoFHandler<dim>::cell_iterator &) const {
+  const vector<vector<Tensor<1,dim> > >&gradients = fe_values.get_shape_grads ();
+  const FullMatrix<double> &values    = fe_values.get_shape_values ();
+  vector<double>            rhs_values (fe_values.n_quadrature_points);
+  const vector<double> &weights   = fe_values.get_JxW_values ();
+  
+  vector<Tensor<1,dim> >   last_solution_grads(fe_values.n_quadrature_points);
+  fe_values.get_function_grads (last_solution, last_solution_grads);
+
+  
+  right_hand_side.value_list (fe_values.get_quadrature_points(), rhs_values);
+   
+  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) += (gradients[i][point] *
+                              gradients[j][point]) *
+                             weights[point] /
+                             sqrt(1+last_solution_grads[i]*last_solution_grads[i]);
+       rhs(i) += values(i,point) *
+                 rhs_values[point] *
+                 weights[point];
+      };
+};
+
+
+
+
+template <int dim>
+void PoissonEquation<dim>::assemble (FullMatrix<double>  &,
+                                    const FEValues<dim> &,
+                                    const DoFHandler<dim>::cell_iterator &) const {
+  Assert (false, ExcPureVirtualFunctionCalled());
+};
+
+
+
+template <int dim>
+void PoissonEquation<dim>::assemble (Vector<double>      &,
+                                    const FEValues<dim> &,
+                                    const DoFHandler<dim>::cell_iterator &) const {
+  Assert (false, ExcPureVirtualFunctionCalled());
+};
+
+
+
+
+template <int dim>
+NonlinearProblem<dim>::NonlinearProblem () :
+               tria(0), dof(0) {};
+
+
+
+template <int dim>
+void NonlinearProblem<dim>::run () {
+
+                                  // first reset everything to a virgin state
+  clear ();
+  
+  tria = new Triangulation<dim>();
+  dof = new DoFHandler<dim> (tria);
+  set_tria_and_dof (tria, dof);
+
+
+  RightHandSide<dim>    rhs;
+  ZeroFunction<dim>     boundary_values;
+  StraightBoundary<dim> boundary;
+  
+  FELinear<dim>                   fe;
+  PoissonEquation<dim>            equation (rhs, last_solution);
+  QGauss2<dim>                    quadrature;
+  
+  ProblemBase<dim>::FunctionMap dirichlet_bc;
+  dirichlet_bc[0] = &boundary_values;
+
+
+  cout << "    Making grid... ";
+  tria->create_hypercube ();
+  tria->refine_global (6);
+  cout << tria->n_active_cells() << " active cells." << endl;
+  
+  cout << "    Distributing dofs... "; 
+  dof->distribute_dofs (fe);
+  cout << dof->n_dofs() << " degrees of freedom." << endl;
+
+                                  // set the starting values for the iteration
+                                  // to a constant value of 1
+  last_solution.reinit (dof->n_dofs());
+  for (unsigned int i=0; i<dof->n_dofs(); ++i)
+    last_solution(i) = 1;
+  
+
+                                  // here comes the fixed point iteration
+  for (unsigned int nonlinear_step=0; nonlinear_step<40; ++nonlinear_step)
+    {
+      cout << "    Nonlinear step " << nonlinear_step << endl;
+      cout << "        Assembling matrices..." << endl;
+      assemble (equation, quadrature, fe,
+               UpdateFlags(update_gradients | update_JxW_values |
+                           update_q_points),
+               dirichlet_bc, boundary);
+
+      cout << "        Solving..." << endl;
+      solve ();
+
+      if (nonlinear_step % 4 == 0)
+       {
+         string filename = "nonlinear.";
+         filename += ('0' + (nonlinear_step/4));
+         filename += ".gnuplot";
+         cout << "        Writing to file <" << filename << ">..." << endl;
+         
+         DataOut<dim> out;
+         ofstream gnuplot(filename.c_str());
+         fill_data (out);
+         out.write_gnuplot (gnuplot);
+         gnuplot.close ();
+       };
+
+      last_solution = solution;
+    };
+  
+  delete dof;
+  delete tria;
+  
+  cout << endl;
+};
+
+
+
+
+int main ()
+{
+  NonlinearProblem<2> problem;
+  problem.run ();
+};

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.