--- /dev/null
+# Listing of Parameters
+# ---------------------
+set Compute error = true
+set Error file format = tex
+set Output error tables = true
+
+# Comma separated list of names for the components. This will be used both for
+# error tables in text format and to output the solution to a file. Note that
+# in the case of a vector function the error name which is used to compute the
+# norm (supposing the type of the other components is 'Add') is the first
+# one.
+set Solution names = u1,u2
+
+# Comma separated version of the same thing as above for the latex version of
+# the table.
+set Solution names for latex = u1,u2
+
+# Comma separated list of table names.
+set Table names = error
+set Write error files = false
+
+
+subsection Dirichlet Data
+ # Any constant used inside the function which is not a variable name.
+ set Function constants =
+
+ # Separate vector valued expressions by ';' as ',' is used internally by the
+ # function parser.
+ set Function expression = sin(2*pi*x)*cos(2*pi*y); cos(2*pi*x)*sin(2*pi*y)
+
+ # The name of the variables as they will be used in the function, separated
+ # by ','.
+ set Variable names = x,y,t
+end
+
+
+subsection Domain Parameters
+ set Input mesh file = square
+ set Input mesh format = ucd
+ set Output mesh file = square_out
+ set Path of domain mesh files = mesh/
+
+ # If this is false, then the input mesh file below is ignored and a
+ # hyper-cube is created.
+ set Read domain mesh from file = false
+
+
+ subsection Grid Out Parameters
+ set Format = none
+
+
+ subsection DX
+ # Write all faces, not only boundary
+ set Write all faces = true
+
+ # Write the mesh connectivity as DX grid cells
+ set Write cells = true
+
+ # If cells are written, additionally write their diameter as data for
+ # visualization
+ set Write diameter = false
+
+ # Write faces of cells. These may be boundary faces or all faces between
+ # mesh cells, according to "Write all faces"
+ set Write faces = false
+
+ # Write the volume of each cell as data
+ set Write measure = false
+ end
+
+ subsection Eps
+ # Azimuth of the viw point, that is, the angle in the plane from the
+ # x-axis.
+ set Azimuth = 30
+
+ # Number of points on boundary edges. Increase this beyond 2 to see
+ # curved boundaries.
+ set Boundary points = 2
+
+ # (2D only) Write cell numbers into the centers of cells
+ set Cell number = false
+
+ # Draw lines with user flag set in different color
+ set Color by flag = false
+
+ # Draw different colors according to grid level.
+ set Color by level = false
+
+ # Elevation of the view point above the xy-plane.
+ set Elevation = 30
+
+ # (2D only) if "Cell number" is true, writenumbers in the form
+ # level.number
+ set Level number = false
+
+ # Width of the lines drawn in points
+ set Line width = 0.5
+
+ # Size of the output in points
+ set Size = 300
+
+ # Depending on this parameter, either thewidth or height of the eps is
+ # scaled to "Size"
+ set Size by = width
+
+ # Write numbers for each vertex
+ set Vertex number = false
+ end
+
+ subsection Gnuplot
+ set Boundary points = 2
+ set Cell number = false
+ end
+
+ subsection Msh
+ set Write faces = false
+ set Write lines = false
+ end
+
+ subsection Ucd
+ set Write faces = false
+ set Write lines = false
+ set Write preamble = true
+ end
+
+ subsection XFig
+ set Boundary = true
+ set Boundary points = 0
+ set Boundary style = 0
+ set Boundary width = 3
+ set Fill style = 20
+ set Level color = false
+ set Level depth = true
+ set Line style = 0
+ set Line width = 1
+ end
+
+ end
+
+end
+
+
+subsection Elastic Moduli
+ # Any constant used inside the functions which is not a variable name.
+ set Function constants =
+
+ # Separate different components expressions by ';' as ',' is used internally
+ # by the function parser.
+ set Row 1 = 0.5; 0.5; 0
+ set Row 2 = 0.5; 0
+ set Row 3 = 0.5
+
+ # The name of the variables as they will be used in the function, separated
+ # by ','.
+ set Variable names = x,y,t
+end
+
+
+subsection Exact Solution
+ # Any constant used inside the function which is not a variable name.
+ set Function constants =
+
+ # Separate vector valued expressions by ';' as ',' is used internally by the
+ # function parser.
+ set Function expression = sin(2*pi*x)*cos(2*pi*y); cos(2*pi*x)*sin(2*pi*y)
+
+ # The name of the variables as they will be used in the function, separated
+ # by ','.
+ set Variable names = x,y,t
+end
+
+
+subsection General Parameters
+ # Determines level of screen output
+ set Console Depth = 10
+
+ # Linear Solver Reduction Tolerance
+ set Linear Solver Reducation Tolerance = 1.0E-12
+
+ # Number of Convergence Cycles Past Initial Refinement
+ set Number of Convergence Cycles = 3
+end
+
+
+subsection Loading Parameters
+ # End time for the simulation
+ set End Time for Simulation = 1.0
+
+ # Number of Loading Cycles to split time interval
+ set Number of Loading Steps = 3
+end
+
+
+subsection Neumann Data
+ # Any constant used inside the function which is not a variable name.
+ set Function constants =
+
+ # Separate vector valued expressions by ';' as ',' is used internally by the
+ # function parser.
+ set Function expression = 0; 0
+
+ # The name of the variables as they will be used in the function, separated
+ # by ','.
+ set Variable names = x,y,t
+end
+
+
+subsection Table 0
+ # Evaluate convergence rates and add a column to the table for each computed
+ # norm.
+ set Add convergence rates = true
+
+ # The extra columns to add to the table.
+ set Extra terms = cells,dofs
+
+ # The caption that will go under the table if we write the file in tex
+ # format. The default value for this object is the same name as the table
+ # name.
+ set Latex table caption = error
+
+ # Each component is separated by a semicolon, and each norm by a comma.
+ # Implemented norms are Linfty, L2, W1infty, H1 and Add, which means that
+ # the norm is added to the previous component. Useful for vector valued
+ # functions.
+ set List of error norms to compute = Linfty, L2, H1 ; Linfty, L2, H1
+end
+
+
+subsection Vector Space Parameters
+ # Boundary indicator, followed by semicolomn and a list of components to
+ # which this boundary conditions apply. More boundary indicators can be
+ # separated by semicolumn. 1:0,1,4 ; 2,4:0,2
+ set Dirichlet boundary map = 0,1,2,3:0,1
+
+ # Ordering of the degrees of freedom: none, comp, cuth.
+ set Dof ordering = cuth
+
+ # The finite element space to use. For vector finite elements use the
+ # notation FESystem[FE_Q(2)^2-FE_DGP(1)] (e.g. Navier-Stokes).
+ set Finite element space = FESystem[FE_Q(1)^2]
+
+ # Degree of the mapping. If 0 is used, then a Cartesian mapping is assumed.
+ set Mapping degree = 1
+
+ # Boundary indicators, followed by semicolomn and a list of components to
+ # which this boundary conditions apply. More boundary indicators can be
+ # separated by semicolumn. 1:0,1,4 ; 2,4:0,2
+ set Neumann boundary map =
+
+ # Boundary indicator, followed by semicolomn and a list of components to
+ # which this boundary conditions apply. More boundary indicators can be
+ # separated by semicolumn. 1:0,1,4 ; 2,4:0,2
+ set Other boundary map = 3:0
+
+
+ subsection Grid Parameters
+ set Bottom fraction = .3
+
+ # If this number is greater than zero, the mesh is distorted upon
+ # refinement in order to disrupt its structureness.
+ set Distortion coefficient = 0
+ set Global refinement = 3
+
+ # A number of zero means no limit.
+ set Max number of cells = 0
+
+ # fixed_number: the Top/Bottom threshold fraction of cells are flagged for
+ # refinement/coarsening. fixed_fraction: the cells whose error is
+ # Top/Bottom fraction of the total are refined/coarsened. optmized: try to
+ # reach optimal error distribution, assuming error is divided by 4 upon
+ # refining. global: refine all cells.
+ set Refinement strategy = fixed_number
+ set Top fraction = .3
+ end
+
+end
+
+
+subsection Yield Stress
+ # Any constant used inside the function which is not a variable name.
+ set Function constants =
+
+ # Separate vector valued expressions by ';' as ',' is used internally by the
+ # function parser.
+ set Function expression = 0
+
+ # The name of the variables as they will be used in the function, separated
+ # by ','.
+ set Variable names = x,y,t
+end
+
+
+subsection f - Body Source
+ # Any constant used inside the function which is not a variable name.
+ set Function constants =
+
+ # Separate vector valued expressions by ';' as ',' is used internally by the
+ # function parser.
+ set Function expression = 8*pi*pi*sin(2*pi*x)*cos(2*pi*y); 8*pi*pi*cos(2*pi*x)*sin(2*pi*y)
+
+ # The name of the variables as they will be used in the function, separated
+ # by ','.
+ set Variable names = x,y,t
+end
+
+
--- /dev/null
+# $Id: Makefile.large 18849 2009-05-15 17:59:46Z bangerth $
+
+# The large projects Makefile looks much like the one for small
+# projects. Basically, only the following seven parameters need to be
+# set by you:
+
+application-name = camclay
+
+# The first denotes the dimension for which the program is to be
+# compiled:
+deal_II_dimension = 2
+
+# The second tells us the name of the executable. It is prefixed by
+# `lib/' to designate its destination directory. Note that the program
+# name depends on the dimension, so you can keep copies for the
+# different dimensions around:
+target = lib/$(application-name)-$(deal_II_dimension)d
+
+# The `debug-mode' variable works as in the small projects Makefile:
+debug-mode = on
+
+# And so does the following variable. You will have to set it to
+# something more reasonable, of course.
+clean-up-files =
+
+# Finally, here is a variable which tells the `run' rule which
+# parameters to pass to the executable. Usually, this will be the name
+# of an input file.
+run-parameters = parameter-file.prm
+
+# Now, this is the last variable you need to set, namely the path to
+# the deal.II toplevel directory:
+D = ../../
+
+
+
+#
+#
+# Usually, you will not need to change anything beyond this point.
+#
+#
+# This tells `make' where to find the global settings and rules:
+include $D/common/Make.global_options
+
+
+# First get a list of files belonging to the project. Include files
+# are expected in `include/', while implementation files are expected
+# in `source/'. Object files are placed into `lib/[123]d', using the
+# same base name as the `.cc' file.
+cc-files = $(shell echo source/*.cc)
+o-files = $(cc-files:source/%.cc=lib/$(deal_II_dimension)d/%.$(OBJEXT))
+go-files = $(cc-files:source/%.cc=lib/$(deal_II_dimension)d/%.g.$(OBJEXT))
+h-files = $(wildcard include/*.h)
+lib-h-files = $(shell echo $D/base/include/base/*.h \
+ $D/lac/include/lac/*.h \
+ $D/deal.II/include/*/*.h)
+
+# As before, define a list of libraries. This, of course depends on
+# the dimension in which we are working:
+libs.g = $(lib-deal2-$(deal_II_dimension)d.g) \
+ $(lib-lac.g) \
+ $(lib-base.g)
+libs.o = $(lib-deal2-$(deal_II_dimension)d.o) \
+ $(lib-lac.o) \
+ $(lib-base.o)
+
+
+
+# Now use the information from above to define the set of libraries to
+# link with and the flags to be passed to the compiler:
+ifeq ($(debug-mode),on)
+ libraries = $(go-files) $(libs.g)
+ flags = $(CXXFLAGS.g)
+else
+ libraries = $(o-files) $(libs.o)
+ flags = $(CXXFLAGS.o)
+endif
+
+
+# Then augment the compiler flags by a specification of the dimension
+# for which the program shall be compiled:
+flags += -Ddeal_II_dimension=$(deal_II_dimension)
+
+
+# The following two rules define how to compile C++ files into object
+# files:
+lib/$(deal_II_dimension)d/%.g.$(OBJEXT) :
+ @echo =====$(application-name)=======$(deal_II_dimension)d====debug=====$(MT)== $(<F)
+ @$(CXX) $(flags) -c $< -o $@
+lib/$(deal_II_dimension)d/%.$(OBJEXT) :
+ @echo =====$(application-name)=======$(deal_II_dimension)d====optimized=$(MT)== $(<F)
+ @$(CXX) $(flags) -c $< -o $@
+
+
+
+# Next define how to link the executable
+$(target)$(EXEEXT) : $(libraries) Makefile
+ @echo =====$(application-name)=======$(deal_II_dimension)d==============$(MT)== Linking $(@F)
+ @$(CXX) -o $@ $(libraries) $(LIBS) $(LDFLAGS)
+
+
+
+# Rule how to run the program
+run: $(target)$(EXEEXT)
+ ./$(target)$(EXEEXT) $(run-parameters)
+
+
+# Rule how to clean up. This is split into several different rules to
+# allow for parallel execution of commands:
+clean: clean-lib clean-data
+ -rm -f *~ */*~ */*/*~ lib/Makefile.dep
+
+clean-lib:
+ -rm -f lib/?d/*.$(OBJEXT) lib/?d/*.g.$(OBJEXT) $(target)$(EXEEXT) lib/TAGS
+
+clean-data:
+ -rm -f $(clean-up-files)
+
+
+# Again tell `make' which rules are not meant to produce files:
+.PHONY: clean clean-data clean-lib run
+
+
+
+# Finally produce the list of dependencies. Note that this time, the
+# object files end up in directories of their own, so we have to
+# modify the output a bit. The file with the dependencies is put into
+# `lib/'.
+lib/Makefile.dep: $(cc-files) $(h-files) $(lib-h-files) Makefile
+ @echo =====$(application-name)=======$(deal_II_dimension)d================== Remaking $@
+ @$D/common/scripts/make_dependencies $(INCLUDE) -Blib $(cc-files) \
+ | $(PERL) -p -e 's!^lib/(.*):!lib/$(deal_II_dimension)d/$$1:!g;' \
+ > $@
+
+include lib/Makefile.dep
+
--- /dev/null
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+Introduction about Cam Clay plasticity models.
--- /dev/null
+<h1>Results</h1>
+
+Astonishing results here.
--- /dev/null
+Solving Cam Clay plasticity model, getting the Nobel Prize.
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id: fe_update_flags.h,v 1.31 2005/10/24 04:33:03 guido Exp $
+// Version: $Name: $
+//
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__assemble_flags_h
+#define __deal2__assemble_flags_h
+
+#include <base/config.h>
+/**
+ * The enum type given to the constructors of LocalAssembleBase objects,
+ * telling those objects which data to assemble on each mesh cell.
+ * When the GlobalAssembler calls the local one, it checks for each flag,
+ * and if it finds one, it assemble the corresponding object.
+ *
+ * By default, all flags are off, i.e. no procedure will be called.
+ *
+ * You can select more than one flag by concatenation
+ * using the bitwise or operator|(AssembleFlags,AssembleFlags).
+ */
+enum AssembleFlags
+ {
+ //! No update
+ assemble_default = 0,
+ //! Cell term.
+ /**
+ * Assemble the cell term. This is
+ * usually needed, unless the matrix
+ * is only a flux matrix.
+ */
+ assemble_cell = 0x0001,
+ //! Assemble boundary term.
+ /**
+ * Calls the assemble_boundary_term for
+ * each boundary face.
+ */
+ assemble_boundary = 0x0002,
+ //! Assemble face term.
+ /** Call the assemble_face_term for each face of each cell in the
+ * triangulation
+ */
+ assemble_face = 0x0004,
+ /** Assemble rhs cell term. Used in assemble_rhs method.*/
+ assemble_rhs_cell = 0x0008,
+ /** Assemble rhs boundary terms. */
+ assemble_rhs_boundary = 0x0010
+ };
+
+
+
+
+
+/**
+ * Global operator which returns an object in which all bits are set
+ * which are either set in the first or the second argument. This
+ * operator exists since if it did not then the result of the bit-or
+ * <tt>operator |</tt> would be an integer which would in turn trigger
+ * a compiler warning when we tried to assign it to an object of type
+ * AssembleFlags.
+ */
+inline
+AssembleFlags
+operator | (AssembleFlags f1, AssembleFlags f2)
+{
+ return static_cast<AssembleFlags> (
+ static_cast<unsigned int> (f1) |
+ static_cast<unsigned int> (f2));
+}
+
+
+
+
+/**
+ * Global operator which sets the bits from the second argument also
+ * in the first one.
+ */
+inline
+AssembleFlags &
+operator |= (AssembleFlags &f1, AssembleFlags f2)
+{
+ f1 = f1 | f2;
+ return f1;
+}
+
+
+/**
+ * Global operator which returns an object in which all bits are set
+ * which are set in the first as well as the second argument. This
+ * operator exists since if it did not then the result of the bit-and
+ * <tt>operator &</tt> would be an integer which would in turn trigger
+ * a compiler warning when we tried to assign it to an object of type
+ * AssembleFlags.
+ */
+inline
+AssembleFlags
+operator & (AssembleFlags f1, AssembleFlags f2)
+{
+ return static_cast<AssembleFlags> (
+ static_cast<unsigned int> (f1) &
+ static_cast<unsigned int> (f2));
+}
+
+
+/**
+ * Global operator which clears all the bits in the first argument if
+ * they are not also set in the second argument.
+ */
+inline
+AssembleFlags &
+operator &= (AssembleFlags &f1, AssembleFlags f2)
+{
+ f1 = f1 & f2;
+ return f1;
+}
+#endif
--- /dev/null
+#ifndef BASE_H
+#define BASE_H
+
+#include <base/function.h>
+#include <base/logstream.h>
+#include <base/utilities.h>
+#include <base/parsed_function.h>
+#include <base/timer.h>
+#include <grid/grid_tools.h>
+#include <grid/grid_generator.h>
+#include <numerics/data_out.h>
+
+#include <grid/filtered_iterator.h>
+
+#include "vector_space.h"
+#include "domain.h"
+#include "error_handler.h"
+#include "linear_elastic.h"
+#include "camclay.h"
+
+
+
+using namespace dealii;
+using namespace dealii::Functions;
+
+template <int dim>
+class Base
+{
+ public:
+ Base ();
+ ~Base ();
+
+ //the routine called from int main()
+ void run ();
+
+ private:
+
+ //function with actual work, called from run()
+ void run_cc (const unsigned int &cc);
+
+ //function with actual work, called from run()
+ void run_step (const unsigned int &cc,
+ const unsigned int &step);
+
+ //parse the parameters and get things rollin
+ void parameters ();
+
+ //elastic predictor
+ void elastic_predictor (double &time);
+
+ //plastic corrector
+ void plastic_corrector ();
+
+ //writes the files
+ void write_files (const unsigned int &cc,
+ const unsigned int &step,
+ const Vector<double> &my_elas,
+ const Vector<double> &my_plas,
+ const Vector<double> &my_other);
+
+ void write_stress_strain (const unsigned int &cc,
+ const unsigned int &step,
+ const Vector<double> &my_elas);
+
+ void write_plot_values(const unsigned int &cc,
+ const unsigned int &step,
+ const Vector<double> &my_elas,
+ const Vector<double> &my_hardening,
+ const Vector<double> &my_plastic,
+ const ParsedSymmetricTensorFunction<4,dim> &C);
+
+ //A parameter handler object
+ ParameterHandler prm;
+
+ //ErrorHandler Object
+ ErrorHandler<dim, Vector<double> > error_handler;
+
+ //This holds the triangulation information
+ Domain<dim> domain;
+
+ //A vectorspace for the thermoelasticity
+ VectorSpace<dim> vspace;
+
+ //Elasticity Class - to change models, just change the class
+ LinearElastic<dim> elastic;
+ //HypoElastic<dim> elastic;
+ //HyperElastic<dim> elastic;
+ //LinearThermoElastic<dim> elastic;
+
+ //Plasticity Class - to change models, just change the class name
+ CamClay<dim> plastic;
+ //CamClayExplicit<dim> plastic;
+ //CamClayInvariant<dim> plastic;
+ //CamClayShearBand<dim> plastic;
+ //J2Flow<dim> plastic;
+ //J2FlowThermal<dim> plastic;
+
+ //General Parameters
+ double lin_red_tol;
+ unsigned int num_steps;
+ double end_time;
+ unsigned int num_cc;
+ unsigned int num_threads;
+ unsigned int console_depth;
+
+ //various functions for the laplace equation
+ ParsedFunction<dim> exact_solution;
+
+
+};
+
+
+#endif
--- /dev/null
+#include "../include/base.h"
+
+template <int dim>
+Base<dim>::Base ()
+ :
+ exact_solution(dim)
+{}
+
+template <int dim>
+Base<dim>::~Base ()
+{}
+
+template <int dim>
+void Base<dim>::run()
+{
+ //setup the depth of screen reporting
+ deallog.depth_console (10);
+
+ //get the parameters
+ parameters ();
+
+ //run as many convergence cycles as we need
+ for(unsigned int cc=0; cc<num_cc; ++cc) {
+ run_cc (cc);
+ }
+
+ //at the end of the problem, output the error table
+ error_handler.output_table(0);
+
+}
+
+template <int dim>
+void Base<dim>::run_cc(const unsigned int &cc)
+{
+ deallog.push("CC-"+Utilities::int_to_string(cc));
+ deallog << "############### Convergence Cycle: "
+ << Utilities::int_to_string(cc)
+ << " ###############" << std::endl;
+
+ //for each convergence cycle, refine
+ //for now, we just leave things global
+ if(cc > 0) vspace.get_tria().refine_global(1);
+
+ vspace.redistribute_dofs();
+
+ elastic.reinit(vspace);
+
+ plastic.reinit(vspace);
+
+ plastic.initial_conditions(vspace);
+
+ deallog << "Building Matrix..." << std::endl;
+ elastic.build_matrix(prm, vspace);
+ plastic.build_matrix(vspace);
+
+ write_files (cc, 0,
+ elastic.sol_total,
+ plastic.solution,
+ plastic.sol_hard_iter);
+
+ write_plot_values(cc, 0, elastic.sol_total, plastic.sol_hard_iter,
+ plastic.solution, elastic.C);
+
+ //just loop over the loading steps
+ for(unsigned int step=1; step<num_steps; ++step) {
+
+ run_step(cc, step);
+
+ }
+
+ //we are interested in the errors, after the whole run
+ exact_solution.set_time(end_time);
+
+ error_handler.error_from_exact(vspace.get_dh(),
+ elastic.sol_total,
+ exact_solution,
+ 0, 0);
+
+ deallog.pop();
+
+}
+
+template <int dim>
+void Base<dim>::run_step(const unsigned int &cc,
+ const unsigned int &step)
+{
+ deallog << std::endl;
+ deallog.push("Step-" + Utilities::int_to_string(step,2));
+
+ double time = end_time * (double(step)/double(num_steps));
+
+ deallog << "Time: " << time <<std::endl;
+
+ elastic_predictor (time);
+
+ plastic_corrector ();
+
+ write_files (cc, step,
+ elastic.sol_total,
+ plastic.solution,
+ plastic.sol_hard_iter);
+
+ write_stress_strain (cc, step, elastic.sol_total);
+
+ write_plot_values(cc, step, elastic.sol_total, plastic.sol_hard_iter,
+ plastic.solution, elastic.C);
+
+ deallog.pop();
+}
+
+template <int dim>
+void Base<dim>::elastic_predictor (double &time)
+{
+ deallog.push("ELASTIC");
+
+ elastic.reinit_step(time);
+
+ elastic.build_rhs(vspace);
+
+ elastic.solve(vspace, lin_red_tol);
+
+ deallog.pop();
+}
+
+template <int dim>
+void Base<dim>::plastic_corrector ()
+{
+ deallog.push("PLASTIC");
+
+ plastic.reinit_step();
+
+ plastic.update_internal_variables(vspace,
+ elastic.sol_total,
+ elastic.C);
+
+ Timer timer;
+ timer.start();
+ plastic.project_strain(prm, vspace, elastic.A);
+ plastic.project_hardening(vspace);
+ timer.stop();
+ deallog << "Time for plastic projection: " << timer() << std::endl;
+
+ deallog.pop();
+}
+
+
+template<int dim>
+void Base<dim>::write_files (const unsigned int &cc,
+ const unsigned int &step,
+ const Vector<double> &my_elas,
+ const Vector<double> &my_plas,
+ const Vector<double> &my_other)
+{
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler (vspace.get_dh());
+
+ std::vector<std::string> elastic_names;
+ std::vector<std::string> plastic_names;
+ std::vector<std::string> other_names;
+ if(dim == 1) {
+ elastic_names.push_back("u1_elastic");
+ plastic_names.push_back("u1_plastic");
+ }
+
+ if(dim == 2) {
+ elastic_names.push_back("u1_elastic");
+ elastic_names.push_back("u2_elastic");
+ plastic_names.push_back("u1_plastic");
+ plastic_names.push_back("u2_plastic");
+ other_names.push_back("hardening");
+ other_names.push_back("Newton_Iterations");
+ }
+
+ if(dim == 3) {
+ elastic_names.push_back("u1_elastic");
+ elastic_names.push_back("u2_elastic");
+ elastic_names.push_back("u3_elastic");
+ plastic_names.push_back("u1_plastic");
+ plastic_names.push_back("u2_plastic");
+ plastic_names.push_back("u3_plastic");
+ other_names.push_back("hardening");
+ other_names.push_back("Newton_Iterations");
+ other_names.push_back("Nothing");
+ }
+
+ data_out.add_data_vector(my_elas, elastic_names);
+ data_out.add_data_vector(my_plas, plastic_names);
+ data_out.add_data_vector(my_other, other_names);
+ data_out.add_data_vector(plastic.iterations, "Average_Newton_Iterations");
+
+ std::vector<unsigned int> partition_int (vspace.get_tria().n_active_cells());
+ GridTools::get_subdomain_association (vspace.get_tria(), partition_int);
+ const Vector<double> partitioning(partition_int.begin(),
+ partition_int.end());
+ data_out.add_data_vector (partitioning, "partitioning");
+
+ data_out.build_patches (vspace.get_fe().degree);
+
+ std::ostringstream filename;
+ filename << "out/solution-";
+ //first the convergence cycle
+ filename << std::setfill('0');
+ filename << std::setw(2) << cc;
+ filename << "-";
+ //then the step
+ filename << std::setfill('0');
+ filename << std::setw(5) << step;
+ filename << ".vtk";
+
+ std::ofstream output (filename.str().c_str());
+ data_out.write_vtk (output);
+
+}
+
+template<int dim>
+void Base<dim>::write_stress_strain (const unsigned int &cc,
+ const unsigned int &step,
+ const Vector<double> &my_elas)
+{
+ //reopen the just written file
+ std::ostringstream filename;
+ filename << "out/solution-";
+ //first the convergence cycle
+ filename << std::setfill('0');
+ filename << std::setw(2) << cc;
+ filename << "-";
+ //then the step
+ filename << std::setfill('0');
+ filename << std::setw(5) << step;
+ filename << ".vtk";
+ std::ofstream output (filename.str().c_str(),std::ios::app);
+
+ //WARNING - THE FOLLOWING IS NOT DIM INDEPENDENT!!!
+
+ //create an fe object to get the gradient values
+ QIterated<dim> qf(QTrapez<1>(),2);
+ FEValues<dim> fe_v(vspace.get_fe(), qf,
+ update_gradients | update_quadrature_points);
+ std::vector< std::vector< Tensor<1,dim> > > total_grads(qf.n_quadrature_points,
+ std::vector<Tensor<1,dim> >(dim));
+ std::vector< Point<dim> > points(qf.n_quadrature_points);
+
+ output << "SCALARS uyy double 1" << std::endl;
+ output << "LOOKUP_TABLE default" << std::endl;
+
+ typename MGDoFHandler<dim>::active_cell_iterator cell = vspace.get_dh().begin_active(),
+ endc = vspace.get_dh().end();
+ for (; cell!=endc; ++cell) {
+
+ fe_v.reinit(cell);
+ fe_v.get_function_gradients(my_elas, total_grads);
+
+ points = fe_v.get_quadrature_points();
+
+ //loop over the quadrature points---------------------------------
+ for (unsigned int qp = 0; qp<qf.n_quadrature_points; ++qp) {
+
+ output << total_grads[qp][1][1] << " ";
+
+ }
+ }
+
+
+ //be responsible - close your files
+ output.close();
+
+}
+
+template<int dim>
+void Base<dim>::write_plot_values(const unsigned int &cc,
+ const unsigned int &step,
+ const Vector<double> &my_elas,
+ const Vector<double> &my_hardening,
+ const Vector<double> &my_plas,
+ const ParsedSymmetricTensorFunction<4,dim> &C)
+{
+ double p=0;
+ double q=0;
+ double deve=0;
+ double k = 0;
+ //std::vector< std::vector< double> > s(dim, std::vector<double>(dim));
+
+ //reopen the just written file
+ std::ostringstream filename;
+ filename << "plot_data/data-";
+ //first the convergence cycle
+ filename << std::setfill('0');
+ filename << std::setw(2) << cc;
+ filename << ".txt";
+ std::ofstream output;
+
+ if(step == 0) {
+ output.open(filename.str().c_str(),std::ios::trunc);
+ } else {
+ output.open(filename.str().c_str(),std::ios::app);
+ }
+
+ //create an fe object to get the gradient values at the nodes
+ //QIterated<dim> qf(QTrapez<1>(),1);
+ QTrapez<dim> qf;
+ FEValues<dim> fe_v(vspace.get_fe(), qf,
+ update_values | update_gradients | update_quadrature_points);
+ std::vector< std::vector< Tensor<1,dim> > > elastic_grads(qf.n_quadrature_points,
+ std::vector<Tensor<1,dim> >(dim));
+ std::vector< std::vector< Tensor<1,dim> > > plastic_grads(qf.n_quadrature_points,
+ std::vector<Tensor<1,dim> >(dim));
+ std::vector< Vector<double> > hard_values(qf.n_quadrature_points, Vector<double>(dim));
+ std::vector< Point<dim> > points(qf.n_quadrature_points);
+
+ //std::vector< std::vector<double> > strain(dim, std::vector<double>(dim) );
+
+ //calculate p,q,deve
+ typename MGDoFHandler<dim>::active_cell_iterator cell = vspace.get_dh().begin_active(),
+ endc = vspace.get_dh().end();
+ for (; cell!=endc; ++cell) {
+
+ fe_v.reinit(cell);
+
+ if (cell->index() == 20 ) {
+ fe_v.get_function_gradients(my_elas, elastic_grads);
+ fe_v.get_function_gradients(my_plas, plastic_grads);
+ fe_v.get_function_values(my_hardening, hard_values);
+
+ points = fe_v.get_quadrature_points();
+
+ SymmetricTensor<2,dim> strain;
+ SymmetricTensor<2,dim> stress;
+ SymmetricTensor<2,3> s;
+
+ for(unsigned int a=0; a<dim; ++a) {
+ for(unsigned int b=0; b<dim; ++b) {
+ strain[a][b] = ( elastic_grads[0][a][b] - plastic_grads[0][a][b]);
+ }
+ }
+
+ for(unsigned int a=0; a<dim; ++a) {
+ for(unsigned int b=0; b<dim; ++b) {
+ for(unsigned int m=0; m<dim; ++m) {
+ for(unsigned int n=0; n<dim; ++n) {
+ stress[a][b] = C(points[0])[a][b][m][n] * strain[m][n];
+ }
+ }
+ }
+ }
+ //explicitly assuming the loading/symmetry condition
+
+ p = (-1.0/3.0) * (2*stress[0][0] + stress[1][1]);
+
+ s[0][0] = stress[0][0];
+ s[1][1] = stress[0][0];
+ s[2][2] = stress[1][1];
+
+ for(unsigned int kk=0; kk<3; ++kk) s[kk][kk] += p;
+
+ q = std::sqrt(3.0/2.0) * s.norm();
+
+ k = hard_values[0](0);
+
+ }
+ }
+
+ //write what we want
+ output << step << " "
+ << p << " "
+ << q << " "
+ << deve << " "
+ << k << " " << std::endl;
+
+
+ //be responsible - close your files
+ output.close();
+
+}
+
+template<int dim>
+void Base<dim>::parameters ()
+{
+ deallog.push("PARAMETERS");
+
+ domain.declare_parameters(prm);
+
+ //add settings for the elastic vector space to the prm file
+ std::string space_name = "Vector Space Parameters";
+ vspace.declare_parameters(prm, space_name);
+
+ error_handler.declare_parameters(prm);
+
+ elastic.declare_parameters(prm);
+
+ plastic.declare_parameters(prm);
+
+ prm.enter_subsection("General Parameters");
+ prm.declare_entry ("Linear Solver Reducation Tolerance",
+ "1.0E-16",
+ Patterns::Double(),
+ "Linear Solver Reduction Tolerance");
+ prm.declare_entry ("Number of Convergence Cycles",
+ "1",
+ Patterns::Integer(),
+ "Number of Convergence Cycles Past Initial Refinement");
+ prm.declare_entry ("Console Depth",
+ "10",
+ Patterns::Integer(),
+ "Determines level of screen output");
+ prm.leave_subsection();
+
+ prm.enter_subsection("Loading Parameters");
+ prm.declare_entry ("End Time for Simulation",
+ "1.0",
+ Patterns::Double(),
+ "End time for the simulation");
+ prm.declare_entry ("Number of Loading Steps",
+ "10",
+ Patterns::Integer(),
+ "Number of Loading Cycles to split time interval");
+ prm.leave_subsection();
+
+ prm.enter_subsection("Exact Solution");
+ Functions::ParsedFunction<dim>::declare_parameters(prm, dim);
+ prm.leave_subsection();
+
+ // ==============================
+
+ prm.read_input("CamClay.prm");
+
+ // ==============================
+ //Initialize the domain
+ domain.reinit(prm);
+ /*
+ domain.get_tria().clear();
+
+ Point<dim> lowerleft(0.0,0.0);
+ Point<dim> upperright(0.3,0.1);
+ std::vector<unsigned int> subs(2);
+ subs[0] = 30; subs[1]=10;
+ GridGenerator::subdivided_hyper_rectangle(domain.get_tria(),
+ subs,
+ lowerleft,
+ upperright, true);
+ */
+ vspace.reinit(prm, domain.get_tria());
+
+ error_handler.parse_parameters(prm);
+
+ elastic.parse_parameters(prm);
+
+ plastic.parse_parameters(prm);
+
+ prm.enter_subsection("Exact Solution");
+ exact_solution.parse_parameters(prm);
+ prm.leave_subsection();
+
+ prm.enter_subsection("General Parameters");
+ lin_red_tol = prm.get_double("Linear Solver Reducation Tolerance");
+ num_cc = prm.get_integer("Number of Convergence Cycles");
+ console_depth = prm.get_integer("Console Depth");
+ prm.leave_subsection();
+
+ prm.enter_subsection("Loading Parameters");
+ num_steps = prm.get_integer("Number of Loading Steps");
+ end_time = prm.get_double("End Time for Simulation");
+ prm.leave_subsection();
+
+ deallog.pop();
+
+}
+
+
--- /dev/null
+#ifndef CAMCLAY
+#define CAMCLAY
+
+//deal.ii packages
+#include <base/logstream.h>
+#include <base/parameter_handler.h>
+#include <base/parsed_function.h>
+#include <base/symmetric_tensor.h>
+#include <base/table.h>
+#include <base/quadrature_lib.h>
+#include <base/utilities.h>
+#include <lac/sparse_direct.h>
+#include <lac/sparse_matrix.h>
+#include <lac/full_matrix.h>
+#include <lac/solver_control.h>
+#include <lac/solver_gmres.h>
+#include <fe/fe.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_tools.h>
+#include <numerics/matrices.h>
+
+//Trilinos packages
+#include <Sacado.hpp>
+
+//LPCM includes
+#include "vector_space.h"
+#include "local_assemble_plastic_project.h"
+#include "local_assemble_scalar_project.h"
+#include "my_tools.h"
+#include "domain.h"
+
+using namespace dealii;
+using namespace dealii::Functions;
+
+template <int dim>
+class CamClay
+{
+ public:
+
+ CamClay();
+
+ ~CamClay();
+
+ void declare_parameters(ParameterHandler &prm);
+
+ void parse_parameters(ParameterHandler &prm);
+
+ void reinit(VectorSpace<dim> &vspace);
+
+ void reinit_step();
+
+ void project_strain(ParameterHandler &prm,
+ VectorSpace<dim> &vspace,
+ SparseMatrix<double> &A);
+
+ void project_hardening(VectorSpace<dim> &vspace);
+
+ void build_matrix(VectorSpace<dim> &vspace);
+
+ void initial_conditions(VectorSpace<dim> &vspace);
+
+ void update_internal_variables(VectorSpace<dim> &vspace,
+ Vector<double> &elastic_solution,
+ ParsedSymmetricTensorFunction<4,dim> &C);
+
+ void compute_jacobian(FullMatrix<double> &jac,
+ const Vector<double> &sol,
+ const double p,
+ const double q,
+ const double M,
+ const SymmetricTensor<4,dim> &C_qp,
+ const SymmetricTensor<2,dim> &xi);
+
+ Vector<double> dF_dstress(const SymmetricTensor<2,dim> &stress,
+ const double k);
+
+ //this is the plastic strain projected into the displacement
+ Vector<double> solution;
+
+ //this is the average number of iterations per cell for the current step
+ Vector<double> iterations;
+
+ //internal variable vectors
+ Table<2, SymmetricTensor<2,dim> > plastic_strain;
+ Table<2, double > hardening;
+ Table<2, double> iter_table;
+
+ Vector<double> sol_hard_iter;
+
+ private:
+
+ //computes the yeild function
+ inline double yield_function(const double p, const double q,
+ const double k);
+
+ //does the actual solving for the new internal variables
+ double solve(const int &index,
+ const unsigned int &qp,
+ const SymmetricTensor<2,dim> &trial_strain);
+
+ //the potential for the hardening
+ double h(const double p,
+ const double k);
+
+ //computes the hydrostatic stress invariant
+ inline double p(const SymmetricTensor<2,dim> &stress);
+
+ //computes the deviatoric stress invariant
+ inline double q(const SymmetricTensor<2,3> &xi);
+
+ //computes the deviatoric stress tensor
+ inline SymmetricTensor<2,3> xi(const SymmetricTensor<2,dim> &stress);
+
+ //computer_theresidua and the jacobian
+ void compute_res_jac(Vector<double> &res,
+ FullMatrix<double> &jac,
+ const Vector<double> &sol,
+ const Vector<double> &prev_step,
+ const SymmetricTensor<2,dim> &total_strain);
+ /*
+ //computes the residual
+ void compute_residual(Vector<double> &res,
+ const Vector<double> &prev_step,
+ const Vector<double> &sol,
+ const SymmetricTensor<2,dim> &stress);
+
+
+
+ //creates the jacobian numerical
+ void numerical_jacobian(FullMatrix<double> &jac,
+ const Vector<double> &sol,
+ const Vector<double> &prev_sol,
+ const SymmetricTensor<2,dim> &total_strain);*/
+
+ //translates from sym tensor notation to voigt notation
+ inline unsigned int sym2voigt(const unsigned int i,
+ const unsigned int j);
+
+ inline unsigned int sac_num(const unsigned int i,
+ const unsigned int j);
+
+ //translates voigt notation to sym tensor indices
+ inline std::vector<unsigned int> voigt2sym(const unsigned int i);
+
+ //delta tensor
+ inline double delta(const unsigned int i,
+ const unsigned int j);
+
+ //voigt delta tensor
+ inline double delta(const unsigned int i);
+
+ //the yeild stress function - not needed for cam clay actually
+ ParsedFunction<dim> yield_stress;
+
+ //some constants
+ //double shear_mod;
+ //double bulk_mod;
+ double M;
+
+ SymmetricTensor<4,dim> C_qp;
+
+ SparseMatrix<double> MM;
+
+ SparsityPattern sp_MM;
+
+
+
+};
+
+#endif
--- /dev/null
+#include "../include/camclay.h"
+
+template <int dim>
+CamClay<dim>::CamClay()
+ :
+yield_stress(1)
+{}
+
+template <int dim>
+CamClay<dim>::~CamClay()
+{
+ MM.clear();
+}
+
+template <int dim>
+void CamClay<dim>::declare_parameters(ParameterHandler &prm)
+{
+ prm.enter_subsection("Yield Stress");
+ Functions::ParsedFunction<dim>::declare_parameters(prm, 1);
+ prm.leave_subsection();
+
+}
+
+template <int dim>
+void CamClay<dim>::parse_parameters(ParameterHandler &prm)
+{
+ prm.enter_subsection("Yield Stress");
+ yield_stress.parse_parameters(prm);
+ prm.leave_subsection();
+}
+
+template <int dim>
+void CamClay<dim>::reinit(VectorSpace<dim> &vspace)
+{
+
+ //the total number of quadrature points
+ double total_qp = std::pow(double(2 * vspace.get_fe().degree + 1), dim);
+
+ //Resize the tables to hold the plastic strains and other internal variables.
+ plastic_strain.reinit(vspace.get_tria().n_active_cells(), int(total_qp));
+ hardening.reinit(vspace.get_tria().n_active_cells(), int(total_qp));
+ iter_table.reinit(vspace.get_tria().n_active_cells(), int(total_qp));
+
+ solution.reinit(vspace.n_dofs());
+ sol_hard_iter.reinit(vspace.n_dofs());
+
+ iterations.reinit(vspace.get_tria().n_active_cells());
+ iterations = 0;
+
+ //setup for projection mass matrix
+ sp_MM.reinit(vspace.n_dofs(),
+ vspace.n_dofs(),
+ vspace.get_dh().max_couplings_between_dofs());
+
+ DoFTools::make_sparsity_pattern (static_cast<DoFHandler<dim> &> (vspace.get_dh()), sp_MM);
+ sp_MM.compress();
+
+ MM.reinit(sp_MM);
+}
+
+template <int dim>
+void CamClay<dim>::reinit_step()
+{
+
+ solution = 0;
+ iterations = 0;
+ sol_hard_iter = 0;
+
+}
+
+template <int dim>
+void CamClay<dim>::initial_conditions(VectorSpace<dim> &vspace)
+{
+
+ QGauss<dim> qf_v((2*vspace.get_fe().degree) + 1);
+
+ FEValues<dim> fe_v (vspace.get_fe(), qf_v,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ const unsigned int n_qp_v = qf_v.n_quadrature_points;
+
+
+ typename MGDoFHandler<dim>::active_cell_iterator cell = vspace.get_dh().begin_active(),
+ endc = vspace.get_dh().end();
+ for (; cell!=endc; ++cell) {
+ fe_v.reinit(cell);
+
+ for (unsigned int qp = 0; qp<n_qp_v; ++qp) {
+
+ hardening(cell->index(), qp) = 2e8;
+
+ //plastic_strain(cell->index(), qp)[0][1] = 0;
+
+
+ }
+ }
+
+
+
+}
+
+template <int dim>
+void CamClay<dim>::update_internal_variables(VectorSpace<dim> &vspace,
+ Vector<double> &elastic_solution,
+ ParsedSymmetricTensorFunction<4, dim> &C)
+{
+
+ //first thing, we need to loop over the cells
+ //and then the quadrature points
+ QGauss<dim> qf_v((2*vspace.get_fe().degree) + 1);
+
+ FEValues<dim> fe_v (vspace.get_fe(), qf_v,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ const unsigned int n_qp_v = qf_v.n_quadrature_points;
+
+ std::vector<Point<dim> > points(n_qp_v);
+
+ SymmetricTensor<2,dim> trial_strain;
+ SymmetricTensor<2,dim> trial_stress;
+ SymmetricTensor<2,dim> total_strain;
+
+ std::vector< std::vector<Tensor<1,dim> > > total_grads(n_qp_v, std::vector<Tensor<1,dim> >(dim));
+
+ //average number of NR Iterations per cell
+ double nr_iterations_ave = 0;
+
+ //This is the material parameter M, for the yield function/surface
+ //right now, I just consider it constant
+ M = 1.0;
+ unsigned int counter = 0;
+ bool marker = false;
+
+ //loop only over the number of cells on this process
+ //we will thread this later
+ typename MGDoFHandler<dim>::active_cell_iterator cell = vspace.get_dh().begin_active(),
+ endc = vspace.get_dh().end();
+ for (; cell!=endc; ++cell) {
+ //reinit with the correct cell
+ fe_v.reinit(cell);
+ int cell_index = cell->index();
+ nr_iterations_ave = 0;
+ marker = false;
+
+ //get the gradients of the PURELY elastic solution
+ fe_v.get_function_gradients(elastic_solution, total_grads);
+
+ //get the quadrature points for the cell
+ points = fe_v.get_quadrature_points();
+
+ //loop over the quadrature points---------------------------------
+ for (unsigned int qp = 0; qp<n_qp_v; ++qp) {
+
+ //get the elastic moduli at this point
+ C_qp = C(points[qp]);
+
+ //calculate the infinitesimal strain
+ for(unsigned int m=0; m<dim; ++m) {
+ for(unsigned int n=0; n<dim; ++n) {
+ total_strain[m][n] = 0.5*(total_grads[qp][m][n]+total_grads[qp][n][m]);
+ }
+ }
+
+ //get the trial strain and hardening. Here we subtract the existing plastic
+ //strain from the new total strain, giving the trial elastic strain
+ //from this step. The trial hardening is simply the exisiting hardening
+ //parameter
+ trial_strain = (total_strain - plastic_strain(cell_index, qp));
+ //trial_strain *= -1.0;
+
+ //calculate the trial stress. This is of course just the contraction
+ //of the elasticity tensor and the trial strain
+ trial_stress = C_qp * trial_strain;
+
+ //if the yield function is greater than zero, solve
+ //for the new values
+ if ( yield_function(p(trial_stress), q(xi(trial_stress)),
+ hardening(cell_index, qp) ) > 0 )
+ {
+ marker = true;
+ //deallog << "Cell: " << cell_index
+ // << " stress01 value: " << trial_stress[0][1] << std::endl;
+
+ nr_iterations_ave += solve(cell_index,
+ qp,
+ total_strain);
+
+
+ }
+
+ } //quad points
+
+ nr_iterations_ave /= n_qp_v;
+
+ iterations(cell_index) = nr_iterations_ave;
+
+ if (marker == true) ++counter;
+
+ } //cells
+
+ deallog << "Plastic Deformation in: "
+ << counter << " of " << vspace.get_tria().n_active_cells() << " cells." << std::endl;
+
+}
+
+template <int dim>
+inline double CamClay<dim>::p(const SymmetricTensor<2, dim> &stress)
+{
+
+ //expicitly assumes specific 2D symmetry - NO SHEAR
+ double pval = (2*stress[0][0] + stress[1][1]);
+
+ pval *= -1.0/3.0;
+
+ return pval;
+
+}
+
+template <int dim>
+inline SymmetricTensor<2,3> CamClay<dim>::xi(const SymmetricTensor<2,dim> &stress)
+{
+
+ SymmetricTensor<2,3> xival;
+
+ xival[0][0] = stress[0][0];
+ xival[1][1] = stress[0][0];
+ xival[2][2] = stress[1][1];
+
+ for(unsigned int i=0; i<3; ++i) xival[i][i] += p(stress);
+
+ return xival;
+}
+
+template <int dim>
+inline double CamClay<dim>::q(const SymmetricTensor<2, 3> &xi)
+{
+
+ double qval = std::sqrt(3.0/2.0) * xi.norm();
+
+ return qval;
+
+}
+
+template <int dim>
+double CamClay<dim>::solve(const int &index,
+ const unsigned int &qp,
+ const SymmetricTensor<2,dim> &total_strain)
+{
+ /*What we have is a system of non-linear ODE that need to
+ be solved for the proper strains, plastic multiplier, and
+ hardening coefficient. The plastic multiplier is relevant only to
+ this function, so we will only see it here*/
+
+ /*This version is simpler - it is NOT using the invariant based
+ methods of Claudio Tamagnini*/
+
+ //these numbers are the
+ //stress components+hardening+plastic multiplier
+ unsigned int size =( (dim*dim) + dim + 4 )/2;
+
+ //we need vectors and a matrix
+ Vector<double> cell_res(size);
+ Vector<double> cell_sol(size);
+ Vector<double> cell_prev(size);
+ FullMatrix<double> cell_jac(size,size);
+
+ //these will be needed later
+ double first_norm = 0;
+ unsigned int iters = 0;
+
+ //create the guess, we start with the previously converged values of
+ //the plastic strain and hardening, and zero for the plastic mult
+ if (dim == 1) {
+ cell_sol(0) = plastic_strain(index, qp)[0][0]; //e_p
+ cell_sol(1) = hardening(index, qp); //k
+ cell_sol(2) = 0; //delta gamma
+ }
+ if (dim == 2) {
+ cell_sol(0) = plastic_strain(index, qp)[0][0]; //e_p 11
+ cell_sol(1) = plastic_strain(index, qp)[1][1]; //e_p 22
+ cell_sol(2) = plastic_strain(index, qp)[0][1]; //e_p 12/21
+ cell_sol(3) = hardening(index, qp); //k
+ cell_sol(4) = 0; //delta gamma
+ }
+ if (dim == 3) {
+ cell_sol(0) = plastic_strain(index, qp)[0][0]; //e_p xx
+ cell_sol(1) = plastic_strain(index, qp)[1][1]; //e_p yy
+ cell_sol(2) = plastic_strain(index, qp)[2][2]; //e_p zz
+ cell_sol(3) = plastic_strain(index, qp)[1][2]; //e_p yz
+ cell_sol(4) = plastic_strain(index, qp)[0][2]; //e_p xz
+ cell_sol(5) = plastic_strain(index, qp)[0][1]; //e_p xy
+ cell_sol(6) = hardening(index, qp); //k
+ cell_sol(7) = 0; //delta gamma
+ }
+ //since the guess of the solution is the previously converged values
+ cell_prev = cell_sol;
+ first_norm = 0;
+ //loop over a newton raphson scheme.
+ //for right now, max nr iterations is set to 10
+ for(unsigned int n=0; n<10; ++n) {
+ cell_res = 0;
+ cell_jac = 0;
+
+ //compute the residual using the guess
+ compute_res_jac(cell_res, cell_jac,
+ cell_sol, cell_prev, total_strain);
+
+ cell_res *= -1.0;
+
+ if ( (index == 20) && (qp == 0)) {
+ deallog << "Cell Index: " << index
+ << " qp: " << qp << " Iters: " << iters
+ << " Rnorm: " << cell_res.linfty_norm() << std::endl;
+
+ // for (unsigned int i=0; i<size; ++i)
+ //deallog << " Residual Comp: " << cell_res(i) << std::endl;
+ }
+
+
+ //check for convergence, for right now just use 1e-10
+ if(n == 0) first_norm = cell_res.linfty_norm();
+ if((cell_res.linfty_norm()/first_norm) < 1e-13) break;
+
+ //if we didn't converge, we are doing an iteration
+ ++iters;
+
+
+ //deallog << "Solution: " << cell_sol(size-2) << std::endl;
+
+ // std::ofstream out("matrix.txt");
+
+ //cell_jac.print_formatted(out,3,true,0,"0",1,0);
+
+ //invert the 5x5 full matrix
+ cell_jac.gauss_jordan();
+
+ //vmult the inverted matrix and add it to the solution
+ cell_jac.vmult_add(cell_sol, cell_res);
+
+ }
+
+ //put the values into the proper places
+ for(unsigned int i=0; i<dim; ++i) {
+ for(unsigned int j=0; j<dim; ++j) {
+ if (i >= j) plastic_strain(index, qp)[i][j] = cell_sol(sym2voigt(i,j));
+ }
+ }
+ hardening(index,qp) = cell_sol(size-2);
+ iter_table(index,qp) = iters;
+
+ return iters;
+}
+
+template <int dim>
+void CamClay<dim>::compute_res_jac(Vector<double> &res,
+ FullMatrix<double> &jac,
+ const Vector<double> &sol,
+ const Vector<double> &prev_sol,
+ const SymmetricTensor<2,dim> &total_strain)
+{
+ //the critical function, computing the residual and the jacobian
+ //basically, if we can assemble the residual using the sacado
+ //doubles, then we have done all of the work. We just need to be careful about
+ //NOT exploiting the fact that the total strain is symmetric.
+
+ typedef Sacado::Fad::DFad<double> fad_double;
+ unsigned int size = sol.size();
+
+ std::vector< std::vector<fad_double> > estrain(dim, std::vector<fad_double>(dim));
+
+ //vectors of the unknowns and independent vars
+ std::vector<fad_double> x(size); //this is ep11,ep22,ep12,k,dgamma
+ std::vector<fad_double> off_diag(((dim*dim)-dim)/2);
+ std::vector<fad_double> R(size); //this is ep11,ep22,ep12,k,dgamma
+
+ for (unsigned int i=0; i<(2*(size-1))-dim; ++i) {
+ //assign the values of the independent variables
+ if( i<size) {
+ x[i] = sol(i);
+ x[i].diff(i,(2*(size-1))-dim);
+ }
+
+ if(i>=size) {
+ off_diag[i-size] = sol(i-3);
+ off_diag[i-size].diff(i,(2*(size-1))-dim);
+ }
+
+ }
+
+ for(unsigned int i=0; i<dim; ++i) {
+ for (unsigned int j=0; j<dim; ++j) {
+
+ if(i >= j) estrain[i][j] = (total_strain[i][j] - x[sym2voigt(i,j)]);
+
+ if(i < j) estrain[i][j] = (total_strain[i][j] - off_diag[sym2voigt(i,j)]);
+
+ }
+ }
+
+ //now we are set, we have estrain, with the right values,
+ //but with different sacado variables, and useable notation
+ std::vector< std::vector<fad_double> > stress(dim, std::vector<fad_double>(dim));
+ fad_double p;
+ p=0;
+ for(unsigned int i=0; i<dim; ++i) {
+ for(unsigned int j=0; j<dim; ++j) {
+ for (unsigned int m=0; m<dim; ++m) {
+ for (unsigned int n=0; n<dim; ++n) {
+
+ stress[i][j] += C_qp[i][j][m][n] * estrain[m][n];
+
+ }
+ }
+
+ if (i==j) p += stress[i][j];
+
+ }
+ }
+
+ p *= (-1.0/double(dim));
+
+ //deallog << "P value: " << p.val() <<std::endl;
+
+ std::vector< std::vector<fad_double> > s(dim, std::vector<fad_double>(dim));
+ fad_double q;
+ fad_double trs;
+ trs = 0;
+
+ for (unsigned int i=0; i<dim; ++i) {
+ for (unsigned int j=0; j<dim; ++j) {
+
+ s[i][j] = stress[i][j];
+
+ if (i == j) {
+
+ s[i][j] += p;
+ trs += s[i][j];
+
+ }
+
+ q += s[i][j]*s[i][j];
+
+ }
+ }
+
+ q = std::sqrt(3.0*q/2.0);
+ //deallog << " s00 value: " << estrain[0][0].val() << std::endl;
+ //deallog << " s01 value: " << estrain[1][1].val() << std::endl;
+ //deallog << " q value: " << q.val() << std::endl;
+ //deallog << " dg value: " << x[size-1].val() << std::endl;
+
+ for(unsigned int i=0; i<(size-2); ++i) {
+
+ R[i] = x[i] - prev_sol(i) - ( x[size-1]*3.0*s[voigt2sym(i)[0]][voigt2sym(i)[1]]/M/M );
+
+ //deallog << R[i] << std::endl;
+
+ if(i<dim) R[i] -= x[size-1]*( ((x[size-2] - 2.0*p)/double(dim)) + (3.0*trs/M/M/double(dim)) );
+
+ }
+
+ //if (R[0].fastAccessDx(0) - R[1].fastAccessDx(1) > 1e-8) deallog << "Shit is Broke" <<std::endl;
+
+ R[size-2] = x[size-2] - prev_sol(size-2) - (16.0*x[size-2]*x[size-1]*(2.0*p - x[size-2]));
+
+ R[size-1] = (q*q/M/M) + (p*(p-x[size-2]));
+
+ for (unsigned int i=0; i<size; ++i) {
+
+ res(i) = R[i].val();
+ //deallog << " R value: " << R[i] << std::endl;
+
+ for (unsigned int j=0; j<size; ++j) {
+
+ jac(i,j) = R[i].fastAccessDx(j);
+
+
+ }//j
+
+ }//i
+
+
+}
+
+
+template <int dim>
+inline std::vector<unsigned int> CamClay<dim>::voigt2sym(const unsigned int i)
+{
+
+ std::vector<unsigned int> return_indices(2);
+
+ if (dim == 2) {
+ if (i == 0) { return_indices[0]=0; return_indices[1]=0; }
+ if (i == 1) { return_indices[0]=1; return_indices[1]=1; }
+ if (i == 2) { return_indices[0]=0; return_indices[1]=1; }
+ }
+
+ if (dim == 3) {
+ //finish this at some time
+ }
+
+ return return_indices;
+
+}
+
+template <int dim>
+inline double CamClay<dim>::delta(const unsigned int i,
+ const unsigned int j)
+{
+ double delta_val=0;
+
+ if(i == j) delta_val = 1;
+
+ return delta_val;
+
+}
+
+template <int dim>
+inline double CamClay<dim>::delta(const unsigned int i)
+{
+ double delta_val=1;
+
+ if(i >= dim) delta_val = 0;
+
+ return delta_val;
+
+}
+
+
+template <int dim>
+inline unsigned int CamClay<dim>::sym2voigt(const unsigned int i,
+ const unsigned int j)
+{
+ unsigned int voigt_value = 0;
+
+ if (dim == 2) {
+
+ if ( (i == 0) && (j == 0) ) voigt_value = 0;
+ if ( (i == 1) && (j == 1) ) voigt_value = 1;
+ if ( (i == 0) && (j == 1) ) voigt_value = 2;
+ if ( (i == 1) && (j == 0) ) voigt_value = 2;
+
+
+ }
+
+ if (dim == 3) {
+ //needs finished
+
+ }
+
+ return voigt_value;
+
+}
+
+template <int dim>
+double CamClay<dim>::h(const double p,
+ const double k)
+{
+
+ double hval = k * ( 2*p - k );
+
+ return hval;
+
+}
+
+template <int dim>
+Vector<double> CamClay<dim>::dF_dstress(const SymmetricTensor<2, dim> &stress,
+ const double k)
+{
+ SymmetricTensor<2,dim> tmp;
+ Vector<double> tmpv(((dim*dim)+dim)/2);
+
+ double p = (1.0/3.0) * first_invariant(stress);
+
+ tmp = stress;
+ for(unsigned int i=0; i<dim; ++i) tmp[i][i] -= p;
+
+ double norm_xi = tmp.norm();
+
+ double q = (std::sqrt(2.0/3.0)) * norm_xi;
+
+ tmp = stress;
+
+ tmp *= (q/M/M) * std::sqrt(8.0/3.0) * (1.0/norm_xi);
+
+ for(unsigned int i=0; i<dim; ++i )
+ tmp[i][i] += (2.0/3.0)*p - (1.0/3.0)*k - (q*p/M/M)*std::sqrt(8.0/3.0)*(1.0/norm_xi);
+
+ tmpv(0) = tmp[0][0];
+ tmpv(1) = tmp[1][1];
+ tmpv(2) = tmp[0][1];
+
+ return tmpv;
+
+}
+
+template <int dim>
+inline double CamClay<dim>::yield_function(const double p, const double q,
+ const double k)
+{
+ double f_value = p*(p-k) + (q*q/M/M);
+
+ return f_value;
+
+}
+
+
+template <int dim>
+void CamClay<dim>::project_strain(ParameterHandler &prm,
+ VectorSpace<dim> &vspace,
+ SparseMatrix<double> &A)
+{
+
+ LocalAssemblePlasticProject<dim> local_plastic_project;
+ local_plastic_project.reinit(vspace.get_fe(), plastic_strain);
+ local_plastic_project.parameters(prm);
+
+ MyTools::assemble_rhs(vspace.get_dh(), vspace.get_hang(),
+ -1, solution, local_plastic_project);
+
+ SparseDirectUMFPACK direct_solver;
+
+ direct_solver.initialize(A);
+
+ direct_solver.solve(solution);
+
+ vspace.get_hang().distribute(solution);
+
+}
+
+template <int dim>
+void CamClay<dim>::build_matrix(VectorSpace<dim> &vspace)
+{
+ QGauss<dim> quad((2*vspace.get_fe().degree) + 1);
+
+ MatrixCreator::create_mass_matrix(vspace.get_dh(), quad, MM);
+
+}
+
+template <int dim>
+void CamClay<dim>::project_hardening(VectorSpace<dim> &vspace)
+{
+ //we need a mass matrix for a scalar variable over the mesh
+ //so, lets build one from scratch
+
+ LocalAssembleScalarProject<dim> local_ass;
+ local_ass.reinit(vspace.get_fe(), hardening, iter_table);
+
+ MyTools::assemble_rhs(vspace.get_dh(), vspace.get_hang(),
+ -1, sol_hard_iter, local_ass);
+
+ SparseDirectUMFPACK direct_solver;
+
+ direct_solver.initialize(MM);
+
+ direct_solver.solve(sol_hard_iter);
+
+ vspace.get_hang().distribute(sol_hard_iter);
+
+ //project the iterations while we are at it
+ //local_ass.reinit(vspace_cc, iter_table);
+ //MyTools::assemble_rhs(vspace_cc.get_dh(), vspace_cc.get_hang(),
+ // -1, sol_iterations, local_ass);
+
+
+
+}
+
--- /dev/null
+#ifndef DOMAIN_H
+#define DOMAIN_H
+#include <fstream>
+
+#include <grid/grid_refinement.h>
+#include <grid/tria.h>
+#include <grid/grid_out.h>
+#include <base/parameter_handler.h>
+#include <base/path_search.h>
+#include <fe/mapping.h>
+#include <dofs/dof_handler.h>
+
+using namespace dealii;
+
+/**
+ Domain object.
+*/
+template <int dim>
+class Domain : public Subscriptor
+{
+ public:
+ /** Empty constructor. */
+ Domain ();
+
+ /** Full constructor. */
+ Domain (ParameterHandler &prm);
+
+ ~Domain ();
+
+ /** Reinit. */
+ void reinit(ParameterHandler &prm);
+
+ /** Read mesh file name, etc. */
+ void parse_parameters(ParameterHandler &prm);
+
+ /** Generate entries in the given parameter file. */
+ static void declare_parameters(ParameterHandler &prm);
+
+ /** Generate the mesh. In this program, the mesh can be read from an
+ input file generated with gmsh (http://www.geuz.org/gmsh/). */
+ void create_mesh ();
+
+ /** Write the mesh. */
+ void output_mesh(std::ostream &out) const;
+
+
+ /** Write the mesh on the file specified by the parameter handler. */
+ void output_mesh() const;
+
+ /** Reference to the triangulation. */
+ inline Triangulation<dim> & get_tria() {
+ Assert(initialized, ExcNotInitialized());
+ return tria;
+ }
+
+ private:
+ bool initialized;
+
+ bool read_mesh;
+
+ PathSearch search_mesh;
+
+ std::string input_mesh_file_name;
+ std::string input_mesh_format;
+ std::string output_mesh_file_name;
+
+ /** Holds the coarse triangulation. */
+ Triangulation<dim> tria;
+ /** Helper class to output the grid. */
+ GridOut gridout;
+};
+#endif
--- /dev/null
+#include "domain.h"
+#include <base/logstream.h>
+#include <base/utilities.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_in.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_tools.h>
+#include <numerics/vectors.h>
+#include <fe/fe.h>
+
+using namespace std;
+
+template <int dim>
+Domain<dim>::Domain() :
+ search_mesh("MESH", 1)
+{
+ initialized = false;
+}
+
+template <int dim>
+Domain<dim>::Domain(ParameterHandler &prm) :
+ search_mesh("MESH", 1)
+{
+ reinit(prm);
+}
+
+template <int dim>
+Domain<dim>::~Domain()
+{
+ // tria.clear();
+}
+
+
+template <int dim>
+void Domain<dim>::reinit(ParameterHandler &prm)
+{
+ deallog.push("DOMAIN");
+ parse_parameters(prm);
+ deallog << "Generating coarse triangulation." << endl;
+ create_mesh();
+ output_mesh();
+ deallog.pop();
+}
+
+
+template <int dim>
+void Domain<dim>::declare_parameters(ParameterHandler &prm)
+{
+ prm.enter_subsection("Domain Parameters");
+ prm.declare_entry ("Read domain mesh from file", "false", Patterns::Bool(),
+ "If this is false, then the input mesh file below is ignored and a hyper-cube is created.");
+ prm.declare_entry ("Path of domain mesh files", "mesh/", Patterns::Anything());
+ prm.declare_entry ("Input mesh file", "square", Patterns::Anything());
+ prm.declare_entry ("Input mesh format", "ucd",
+ Patterns::Selection(GridIn<dim>::get_format_names()));
+ prm.declare_entry ("Output mesh file", "square_out", Patterns::Anything());
+
+ prm.enter_subsection("Grid Out Parameters");
+ GridOut::declare_parameters(prm);
+ prm.leave_subsection();
+
+ prm.leave_subsection();
+
+}
+
+ template <int dim>
+void Domain<dim>::parse_parameters(ParameterHandler &prm)
+{
+ prm.enter_subsection("Domain Parameters");
+ read_mesh = prm.get_bool ("Read domain mesh from file");
+ search_mesh.add_path(prm.get ("Path of domain mesh files"));
+ input_mesh_file_name = prm.get ("Input mesh file");
+ output_mesh_file_name = prm.get ("Output mesh file");
+ input_mesh_format = prm.get ("Input mesh format");
+
+ prm.enter_subsection("Grid Out Parameters");
+ gridout.parse_parameters(prm);
+ prm.leave_subsection();
+
+ prm.leave_subsection();
+}
+
+template <int dim>
+void Domain<dim>::create_mesh()
+{
+ if(read_mesh) {
+ GridIn<dim> grid_in;
+ grid_in.attach_triangulation (tria);
+ string mfilen = search_mesh.find
+ (input_mesh_file_name,
+ grid_in.default_suffix(grid_in.parse_format(input_mesh_format)),
+ "r");
+ ifstream mfile(mfilen.c_str());
+ grid_in.read(mfile, GridIn<dim>::parse_format(input_mesh_format));
+ } else {
+ Point<dim> corner;
+ for(unsigned int d=0; d<dim; ++d) corner[d] = 1.;
+ GridGenerator::hyper_rectangle (tria, Point<dim>(), corner, true);
+// GridTools::partition_triangulation (n_mpi_processes, tria);
+ }
+ initialized = true;
+}
+
+
+template <int dim>
+void Domain<dim>::output_mesh(std::ostream &out) const {
+ Assert(initialized, ExcNotInitialized());
+ gridout.write (tria, out);
+}
+
+
+template <int dim>
+void Domain<dim>::output_mesh() const {
+ if(gridout.default_suffix() != "") {
+ std::ofstream out_file ((output_mesh_file_name + gridout.default_suffix()).c_str());
+ output_mesh(out_file);
+ out_file.close();
+ }
+}
--- /dev/null
+#ifndef LH_ERROR_HANDLER_H
+#define LH_ERROR_HANDLER_H
+
+#include <fstream>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_constraints.h>
+#include <lac/vector.h>
+
+#include <grid/tria.h>
+
+// #include <numerics/error_estimator.h>
+#include <base/function.h>
+#include <numerics/solution_transfer.h>
+#include <numerics/data_out.h>
+
+#include <base/convergence_table.h>
+#include <base/logstream.h>
+#include <base/config.h>
+
+#include <base/parameter_handler.h>
+
+#include <map>
+
+enum NormFlags {
+ None = 0x00,
+ Linfty = 0x01,
+ L2 = 0x02,
+ W1infty = 0x04,
+ H1 = 0x08,
+ AddUp = 0x10
+};
+
+using namespace dealii;
+
+template <int dim, typename VEC=Vector<double> >
+class ErrorHandler : public Subscriptor
+{
+public:
+ /** The constructor takes the mpi initialization stuff. */
+ ErrorHandler ();
+
+ /** Initialize the given values for the paramter file. */
+ static void declare_parameters(ParameterHandler &prm,
+ unsigned int ntables=1);
+
+ /** Parse the given parameter handler. */
+ void parse_parameters(ParameterHandler &prm);
+
+ /** Calculate the error of the numeric solution in variuous norms. Store
+ the result in the given table. */
+ void error_from_exact(const DoFHandler<dim> & vspace,
+ const VEC &solution,
+ const Function<dim> &exact,
+ unsigned int table_no = 0,
+ double dt=0.);
+
+ /** Difference between two solutions in two different vector spaces. */
+ void difference(const DoFHandler<dim> &, const VEC &,
+ const DoFHandler<dim> &, const VEC &,
+ unsigned int table_no = 0, double dt=0.);
+
+ /** Difference between two solutions in the same vector space. */
+ void difference(const DoFHandler<dim> &, const VEC &,
+ const VEC &, unsigned int table_no = 0, double dt=0.);
+
+ /** By default output first table. */
+ void output_table(const unsigned int table_no=0);
+
+private:
+ /** Error results.*/
+ std::vector<ConvergenceTable> tables;
+
+ /** Headers for tables and output. Contains the name of the solution
+ components. */
+ std::vector<std::string> headers;
+
+ /** Headers for latex tables. Contains the name of the solution
+ components. */
+ std::vector<std::string> latex_headers;
+
+ /** Captions for latex. */
+ std::vector<std::string> latex_captions;
+
+ /** Names of the tables. */
+ std::vector<std::string> names;
+
+ /** Type of error to compute per components. */
+ std::vector<std::vector<NormFlags> > types;
+
+ /** The parameters have been read. */
+ bool initialized;
+
+ /** Compute the error. If this is false, all functions regarding
+ errors are disabled and don't do anything.*/
+ bool compute_error;
+
+ /** Add convergence rates. */
+ std::vector<bool> add_rates;
+
+ /** Write the error files. */
+ bool write_error;
+
+ /** Output the error file also on screen. */
+ bool output_error;
+
+ /** The error file format. */
+ std::string error_file_format;
+
+ /** The extra column to add to the tables. */
+ std::vector<std::map<std::string, bool> > extras;
+
+ /** Wether or not to calculate the rates according to the given keys. */
+ std::vector<std::string> rate_keys;
+};
+
+/**
+ * Global operator which returns an object in which all bits are set
+ * which are either set in the first or the second argument. This
+ * operator exists since if it did not then the result of the bit-or
+ * <tt>operator |</tt> would be an integer which would in turn trigger
+ * a compiler warning when we tried to assign it to an object of type
+ * NormFlags.
+ */
+inline
+NormFlags
+operator | (NormFlags f1, NormFlags f2)
+{
+ return static_cast<NormFlags> (
+ static_cast<unsigned int> (f1) |
+ static_cast<unsigned int> (f2));
+}
+
+/**
+ * Global operator which sets the bits from the second argument also
+ * in the first one.
+ */
+inline
+NormFlags &
+operator |= (NormFlags &f1, NormFlags f2)
+{
+ f1 = f1 | f2;
+ return f1;
+}
+
+
+/**
+ * Global operator which returns an object in which all bits are set
+ * which are set in the first as well as the second argument. This
+ * operator exists since if it did not then the result of the bit-and
+ * <tt>operator &</tt> would be an integer which would in turn trigger
+ * a compiler warning when we tried to assign it to an object of type
+ * NormFlags.
+ */
+inline
+NormFlags
+operator & (NormFlags f1, NormFlags f2)
+{
+ return static_cast<NormFlags> (
+ static_cast<unsigned int> (f1) &
+ static_cast<unsigned int> (f2));
+}
+
+
+/**
+ * Global operator which clears all the bits in the first argument if
+ * they are not also set in the second argument.
+ */
+inline
+NormFlags &
+operator &= (NormFlags &f1, NormFlags f2)
+{
+ f1 = f1 & f2;
+ return f1;
+}
+
+#endif
--- /dev/null
+#include "../include/error_handler.h"
+
+#include <base/logstream.h>
+#include <base/quadrature_lib.h>
+#include <base/utilities.h>
+
+#include <grid/grid_tools.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <fe/mapping_q.h>
+#include <fe/fe.h>
+#include <cstdio>
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <string>
+
+template <int dim, typename VECTOR>
+ErrorHandler<dim,VECTOR>::ErrorHandler ()
+{
+ initialized = false;
+}
+
+template <int dim, typename VECTOR>
+void ErrorHandler<dim,VECTOR>::declare_parameters (ParameterHandler &prm,
+ unsigned int ntables)
+{
+ prm.declare_entry ("Write error files", "false", Patterns::Bool());
+ prm.declare_entry ("Output error tables", "true", Patterns::Bool());
+ prm.declare_entry ("Error file format", "tex", Patterns::Selection("tex|txt"));
+ prm.declare_entry ("Compute error", "true", Patterns::Bool());
+ prm.declare_entry ("Table names", "error", Patterns::Anything(),
+ "Comma separated list of table names. ");
+ prm.declare_entry ("Solution names", "u", Patterns::Anything(),
+ "Comma separated list of names for the components. This "
+ "will be used both for error tables in text format and to "
+ "output the solution to a file. Note that in the case "
+ "of a vector function the error name which is used to "
+ "compute the norm (supposing the type of the other "
+ "components is 'Add') is the first one.");
+ prm.declare_entry ("Solution names for latex", "u", Patterns::Anything(),
+ "Comma separated version of the same thing as above for "
+ "the latex version of the table.");
+
+ // prm.declare_entry ("Ib output format", "msh", Patterns::Selection("raw|msh"));
+ // prm.declare_entry ("Ib input file prefix", "ellipse", Patterns::Anything());
+ for(unsigned int i=0; i<ntables; ++i) {
+ char tmp[10];
+ sprintf(tmp, "Table %d", i);
+ prm.enter_subsection(tmp);
+
+ prm.declare_entry("List of error norms to compute", "Linfty, L2, H1",
+ Patterns::Anything(), "Each component is separated by a semicolon, "
+ "and each norm by a comma. Implemented norms are Linfty, L2, "
+ "H1 and AddUp, which means that the norm is added to the previous "
+ "component. Useful for vector valued functions.");
+ prm.declare_entry("Add convergence rates", "true", Patterns::Bool(),
+ "Evaluate convergence rates and add a column to the table for each "
+ "computed norm. ");
+ prm.declare_entry("Latex table caption", "error", Patterns::Anything(),
+ "The caption that will go under the table if we write the file in "
+ "tex format. The default value for this object is the same name "
+ "as the table name.");
+ prm.declare_entry("Extra terms", "cells,dofs",
+ Patterns::Anything(),
+ "The extra columns to add to the table.");
+ prm.declare_entry("Rate key", "",
+ Patterns::Selection("dofs|cells|dt|"),
+ "The key to use to compute the convergence rates.");
+ prm.leave_subsection();
+ }
+}
+
+template <int dim, typename VECTOR>
+void ErrorHandler<dim,VECTOR>::parse_parameters (ParameterHandler &prm)
+{
+ write_error = prm.get_bool ("Write error files");
+ output_error = prm.get_bool ("Output error tables");
+
+ error_file_format = prm.get ("Error file format");
+ compute_error = prm.get_bool ("Compute error");
+ std::string all_names = prm.get ("Table names");
+ headers = Utilities::split_string_list(prm.get ("Solution names"));
+ latex_headers = Utilities::split_string_list(prm.get ("Solution names for latex"));
+
+ if (all_names != "") {
+ names = Utilities::split_string_list(all_names);
+ types.resize(names.size(), std::vector<NormFlags> (headers.size()));
+ add_rates.resize(names.size());
+ tables.resize(names.size());
+ latex_captions.resize(names.size());
+ std::map<std::string, bool> extra;
+ extra["dof"] = false;
+ extra["cells"] = false;
+ extra["dt"] = false;
+
+ extras.resize(names.size(), extra);
+ rate_keys.resize(names.size(), "");
+
+ for(unsigned int i=0; i<names.size(); ++i) {
+ char tmp[10];
+ sprintf(tmp, "Table %d", i);
+ prm.enter_subsection(tmp);
+
+ all_names = prm.get("List of error norms to compute");
+ add_rates[i] = prm.get_bool("Add convergence rates");
+ rate_keys[i] = prm.get("Rate key");
+ latex_captions[i] = prm.get("Latex table caption");
+ std::vector<std::string> all_extras =
+ Utilities::split_string_list(prm.get("Extra terms"));
+
+ for(unsigned int x=0; x< all_extras.size(); ++x)
+ extras[i][all_extras[x]] = true;
+
+ prm.leave_subsection();
+
+ std::vector<std::string> all_comps = Utilities::split_string_list(all_names, ';');
+ // Check that the input string has all the needed fields
+ AssertThrow(all_comps.size() == headers.size(),
+ ExcDimensionMismatch(all_comps.size() , headers.size()));
+
+ for(unsigned int j=0; j<all_comps.size(); ++j) {
+ std::vector<std::string> all_types =
+ Utilities::split_string_list(all_comps[j]);
+ for(unsigned int k=0; k<all_types.size(); ++k) {
+ if(all_types[k] == "Linfty") {
+ types[i][j] |= Linfty;
+ } else if(all_types[k] == "L2") {
+ types[i][j] |= L2;
+ } else if(all_types[k] == "W1infty") {
+ types[i][j] |= W1infty;
+ } else if(all_types[k] == "H1") {
+ types[i][j] |= H1;
+ } else if(all_types[k] == "AddUp") {
+ types[i][j] |= AddUp;
+ } else {
+ AssertThrow(false, ExcMessage("Didn't recognize a norm type."));
+ }
+ }
+ }
+ }
+ }
+ initialized = true;
+}
+
+template <int dim, typename VECTOR>
+void ErrorHandler<dim,VECTOR>::output_table (const unsigned int table_no) {
+ if (compute_error) {
+ AssertThrow(initialized, ExcNotInitialized());
+ AssertThrow(table_no < names.size(), ExcIndexRange(table_no, 0, names.size()));
+
+ // Add convergence rates
+ if(add_rates[table_no]) {
+ if(extras[table_no]["dofs"])
+ tables[table_no].omit_column_from_convergence_rate_evaluation("dofs");
+ if(extras[table_no]["cells"])
+ tables[table_no].omit_column_from_convergence_rate_evaluation("cells");
+ if(extras[table_no]["dt"])
+ tables[table_no].omit_column_from_convergence_rate_evaluation("dt");
+ if(rate_keys[table_no] == "")
+ tables[table_no].evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2);
+ else
+ tables[table_no].evaluate_all_convergence_rates(rate_keys[table_no], ConvergenceTable::reduction_rate_log2);
+ }
+
+ if(output_error) tables[table_no].write_text(std::cout);
+
+ if(write_error) {
+ std::string filename = names[table_no] +
+ "." + error_file_format;
+
+ std::ofstream table_file(filename.c_str());
+
+ if(error_file_format != "txt")
+ tables[table_no].write_tex(table_file);
+ else
+ tables[table_no].write_text(table_file);
+ table_file.close();
+ }
+ }
+}
+
+template <int dim, typename VECTOR>
+void ErrorHandler<dim,VECTOR>::difference(const DoFHandler<dim> & dh,
+ const VECTOR &solution1,
+ const VECTOR &solution2,
+ unsigned int table_no,
+ double dt) {
+ AssertThrow(solution1.size() == solution2.size(), ExcDimensionMismatch(
+ solution1.size(), solution2.size()));
+ VECTOR solution(solution1);
+ solution -= solution2;
+ error_from_exact(dh, solution,
+ ConstantFunction<dim>(0, headers.size()), table_no, dt);
+}
+
+
+template <int dim, typename VECTOR>
+void ErrorHandler<dim,VECTOR>::error_from_exact(const DoFHandler<dim> & dh,
+ const VECTOR &solution,
+ const Function<dim> &exact,
+ unsigned int table_no,
+ double dt)
+{
+ if (compute_error) {
+ AssertThrow(initialized, ExcNotInitialized());
+ AssertThrow(table_no < types.size(), ExcIndexRange(table_no, 0, names.size()));
+ AssertThrow(exact.n_components == types[table_no].size(),
+ ExcDimensionMismatch(exact.n_components, types[table_no].size()));
+
+ deallog.push("Error");
+ deallog << "Calculating Errors." << std::endl;
+ std::vector< std::vector<double> > error( exact.n_components, std::vector<double>(4));
+ const unsigned int n_active_cells = dh.get_tria().n_active_cells();
+ const unsigned int n_dofs=dh.n_dofs();
+
+ if(extras[table_no]["cells"]) {
+ tables[table_no].add_value("cells", n_active_cells);
+ tables[table_no].set_tex_caption("cells", "\\# cells");
+ tables[table_no].set_tex_format("cells", "r");
+ }
+ if(extras[table_no]["dofs"]) {
+ tables[table_no].add_value("dofs", n_dofs);
+ tables[table_no].set_tex_caption("dofs", "\\# dofs");
+ tables[table_no].set_tex_format("dofs", "r");
+ }
+ if(extras[table_no]["dt"]) {
+ tables[table_no].add_value("dt", dt);
+ tables[table_no].set_tex_caption("dt", "\\Delta t");
+ tables[table_no].set_tex_format("dt", "r");
+ }
+
+ bool compute_Linfty = false;
+ bool compute_L2 = false;
+ bool compute_W1infty = false;
+ bool compute_H1 = false;
+ bool add_this = false;
+
+ unsigned int last_non_add = 0;
+
+ for(unsigned int component=0; component < exact.n_components; ++component) {
+ NormFlags norm = types[table_no][component];
+
+ deallog << "Error flags: " << norm << std::endl;
+
+ // Select one Component
+ ComponentSelectFunction<dim> select_component ( component, 1. , exact.n_components);
+
+ Vector<float> difference_per_cell (dh.get_tria().n_active_cells());
+
+ QGauss<dim> q_gauss((dh.get_fe().degree+1) * 2);
+
+ // The add bit is set
+ add_this = (norm & AddUp);
+
+ if(!add_this) {
+ last_non_add = component;
+ compute_L2 = ( norm & L2 );
+ compute_H1 = ( norm & H1 );
+ compute_W1infty = ( norm & W1infty ) ;
+ compute_Linfty = ( norm & Linfty );
+ }
+ // if add is set, we do not modify the previous selection
+
+ if(compute_L2) {
+ VectorTools::integrate_difference (//mapping,
+ dh, //dof_handler,
+ solution,
+ exact,
+ difference_per_cell,
+ q_gauss,
+ VectorTools::L2_norm,
+ &select_component );
+ }
+
+ const double L2_error = difference_per_cell.l2_norm();
+ difference_per_cell = 0;
+
+ if(compute_H1) {
+ VectorTools::integrate_difference (//mapping,
+ dh, //dof_handler,
+ solution,
+ exact,
+ difference_per_cell,
+ q_gauss,
+ VectorTools::H1_norm,
+ &select_component );
+ }
+ const double H1_error = difference_per_cell.l2_norm();
+ difference_per_cell = 0;
+
+ if(compute_W1infty) {
+ VectorTools::integrate_difference (//mapping,
+ dh, //dof_handler,
+ solution,
+ exact,
+ difference_per_cell,
+ q_gauss,
+ VectorTools::W1infty_norm,
+ &select_component );
+ }
+
+ const double W1inf_error = difference_per_cell.linfty_norm();
+
+ if(compute_Linfty) {
+ VectorTools::integrate_difference (//mapping,
+ dh, //dof_handler,
+ solution,
+ exact,
+ difference_per_cell,
+ q_gauss,
+ VectorTools::Linfty_norm,
+ &select_component );
+ }
+
+ const double Linf_error = difference_per_cell.linfty_norm();
+
+ if(add_this) {
+ AssertThrow(component, ExcMessage("Cannot add on first component!"));
+
+ error[last_non_add][0] = std::max(error[last_non_add][0], Linf_error);
+ error[last_non_add][1] += L2_error;
+ error[last_non_add][2] = std::max(error[last_non_add][2], W1inf_error);
+ error[last_non_add][3] += H1_error;
+
+ } else {
+
+ error[component][0] = Linf_error;
+ error[component][1] = L2_error;
+ error[component][2] = W1inf_error;
+ error[component][3] = H1_error;
+
+ }
+ }
+
+ for(unsigned int j=0; j<exact.n_components; ++j) {
+ NormFlags norm = types[table_no][j];
+ // If this was added, don't do anything
+ if( !(norm & AddUp) ) {
+ if(norm & Linfty) {
+ std::string name = headers[j] + "_Linfty";
+ std::string latex_name = "$\\| " +
+ latex_headers[j] + " - " +
+ latex_headers[j] +"_h \\|_\\infty $";
+ double this_error = error[j][0];
+
+ tables[table_no].add_value(name, this_error);
+ tables[table_no].set_precision(name, 3);
+ tables[table_no].set_scientific(name, true);
+ tables[table_no].set_tex_caption(name, latex_name);
+ }
+
+ if(norm & L2) {
+ std::string name = headers[j] + "_L2";
+ std::string latex_name = "$\\| " +
+ latex_headers[j] + " - " +
+ latex_headers[j] +"_h \\|_0 $";
+ double this_error = error[j][1];
+
+ tables[table_no].add_value(name, this_error);
+ tables[table_no].set_precision(name, 3);
+ tables[table_no].set_scientific(name, true);
+ tables[table_no].set_tex_caption(name, latex_name);
+ }
+ if(norm & W1infty) {
+ std::string name = headers[j] + "_W1infty";
+ std::string latex_name = "$\\| " +
+ latex_headers[j] + " - " +
+ latex_headers[j] +"_h \\|_{1,\\infty} $";
+ double this_error = error[j][2];
+
+ tables[table_no].add_value(name, this_error);
+ tables[table_no].set_precision(name, 3);
+ tables[table_no].set_scientific(name, true);
+ tables[table_no].set_tex_caption(name, latex_name);
+ }
+ if(norm & H1){
+ std::string name = headers[j] + "_H1";
+ std::string latex_name = "$\\| " +
+ latex_headers[j] + " - " +
+ latex_headers[j] +"_h \\|_1 $";
+ double this_error = error[j][3];
+
+ tables[table_no].add_value(name, this_error);
+ tables[table_no].set_precision(name, 3);
+ tables[table_no].set_scientific(name, true);
+ tables[table_no].set_tex_caption(name, latex_name);
+ }
+ }
+ }
+ deallog.pop();
+ }
+}
+
--- /dev/null
+#ifndef LINEARELASTIC
+#define LINEARELASTIC
+
+#include <base/logstream.h>
+#include <base/parameter_handler.h>
+#include <base/parsed_function.h>
+#include <lac/vector.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/sparse_matrix.h>
+#include <lac/filtered_matrix.h>
+#include <lac/precondition.h>
+#include <lac/sparse_ilu.h>
+#include <lac/solver_cg.h>
+#include <lac/sparse_direct.h>
+#include <dofs/dof_tools.h>
+
+
+#include "parsed_symmetric_tensor_function.h"
+#include "local_assemble_elastic_matrix.h"
+#include "local_assemble_elastic_rhs.h"
+#include "vector_space.h"
+#include "my_tools.h"
+
+using namespace dealii;
+using namespace dealii::Functions;
+
+template <int dim>
+class LinearElastic
+{
+ public:
+
+ LinearElastic();
+
+ ~LinearElastic();
+
+ void declare_parameters(ParameterHandler &prm);
+
+ void parse_parameters(ParameterHandler &prm);
+
+ void reinit(VectorSpace<dim> &vspace);
+
+ void build_matrix(ParameterHandler &prm,
+ VectorSpace<dim> &vspace);
+
+ void reinit_step(double &time);
+
+ void build_rhs(VectorSpace<dim> &vspace);
+
+ void solve(VectorSpace<dim> &vspace,
+ double &tolerance);
+
+ SparseMatrix<double> A;
+
+ ParsedSymmetricTensorFunction<4, dim> C;
+
+ //The solution vectors
+ Vector<double> sol_total;
+ Vector<double> sol_increment;
+
+ private:
+
+ //various functions for the laplace equation
+ ParsedFunction<dim> dbc;
+ ParsedFunction<dim> nbc;
+ ParsedFunction<dim> bf;
+
+ //The system matrix for linear elasticity
+ SparsityPattern sp_A;
+
+ //The distributed right hand side
+ Vector<double> rhs;
+
+
+
+};
+
+#endif
--- /dev/null
+#include "../include/linear_elastic.h"
+
+template <int dim>
+LinearElastic<dim>::LinearElastic()
+ :
+ dbc(dim),
+ nbc(dim),
+ bf(dim)
+{}
+
+
+template <int dim>
+LinearElastic<dim>::~LinearElastic()
+{
+ A.clear();
+}
+
+template <int dim>
+void LinearElastic<dim>::declare_parameters(ParameterHandler &prm)
+{
+ prm.enter_subsection("Dirichlet Data");
+ Functions::ParsedFunction<dim>::declare_parameters(prm, dim);
+ prm.leave_subsection();
+
+ prm.enter_subsection("f - Body Source");
+ Functions::ParsedFunction<dim>::declare_parameters(prm, dim);
+ prm.leave_subsection();
+
+ prm.enter_subsection("Elastic Moduli");
+ ParsedSymmetricTensorFunction<4, dim>::declare_parameters(prm);
+ prm.leave_subsection();
+
+ prm.enter_subsection("Neumann Data");
+ Functions::ParsedFunction<dim>::declare_parameters(prm, dim);
+ prm.leave_subsection();
+
+}
+
+template <int dim>
+void LinearElastic<dim>::parse_parameters(ParameterHandler &prm)
+{
+ prm.enter_subsection("f - Body Source");
+ bf.parse_parameters(prm);
+ prm.leave_subsection();
+
+ prm.enter_subsection("Neumann Data");
+ nbc.parse_parameters(prm);
+ prm.leave_subsection();
+
+ prm.enter_subsection("Elastic Moduli");
+ C.parse_parameters(prm);
+ prm.leave_subsection();
+
+ prm.enter_subsection("Dirichlet Data");
+ dbc.parse_parameters(prm);
+ prm.leave_subsection();
+
+}
+
+template <int dim>
+void LinearElastic<dim>::reinit(VectorSpace<dim> &vspace)
+{
+
+ sp_A.reinit(vspace.n_dofs(),
+ vspace.n_dofs(),
+ vspace.get_dh().max_couplings_between_dofs());
+
+ DoFTools::make_sparsity_pattern (static_cast<DoFHandler<dim> &> (vspace.get_dh()),
+ sp_A);
+ vspace.get_hang().condense (sp_A);
+ sp_A.compress();
+
+ A.reinit(sp_A);
+
+ sol_total.reinit(vspace.n_dofs());
+ sol_increment.reinit(vspace.n_dofs());
+ rhs.reinit(vspace.n_dofs());
+
+}
+
+template <int dim>
+void LinearElastic<dim>::build_matrix(ParameterHandler &prm,
+ VectorSpace<dim> &vspace)
+{
+
+ LocalAssembleElasticMatrix<dim> local_elastic_matrix;
+ local_elastic_matrix.reinit(vspace.get_fe());
+ local_elastic_matrix.parameters(prm);
+
+ MyTools::assemble(vspace.get_dh(), vspace.get_hang(),
+ -1, A, local_elastic_matrix);
+
+}
+
+template <int dim>
+void LinearElastic<dim>::reinit_step(double &time)
+{
+
+ rhs = 0;
+ //sol_increment = 0;
+ bf.set_time(time);
+ nbc.set_time(time);
+ dbc.set_time(time);
+}
+
+template <int dim>
+void LinearElastic<dim>::build_rhs(VectorSpace<dim> &vspace)
+{
+
+ LocalAssembleElasticRHS<dim> elastic_local_rhs;
+
+ elastic_local_rhs.reinit(vspace.get_fe(),
+ bf, nbc,
+ vspace.neumann_bc);
+
+ MyTools::assemble_rhs(vspace.get_dh(), vspace.get_hang(),
+ -1, rhs, elastic_local_rhs);
+
+}
+
+template <int dim>
+void LinearElastic<dim>::solve(VectorSpace<dim> &vspace,
+ double &tolerance)
+{
+
+ //get the DBC map
+ std::map<unsigned int, double> bv;
+ vspace.interpolate_dirichlet_bc(dbc, bv);
+
+ //we want to leave the matrix untouched for future projections
+ //so we use a filtered matrix object
+ FilteredMatrix<Vector<double> > filtered_A(A);
+ filtered_A.add_constraints(bv);
+ filtered_A.apply_constraints(rhs, true);
+
+ deallog << "RHS L2 Norm After DBC: " << rhs.l2_norm() << std::endl;
+
+ //make the preconditioner
+ PreconditionJacobi<SparseMatrix<double> > precon;
+ precon.initialize(A, 0.8);
+ FilteredMatrix<Vector<double> > filtered_precon(precon);
+
+ SolverControl control (vspace.n_dofs(),
+ tolerance*rhs.linfty_norm(),
+ false, true);
+
+ GrowingVectorMemory<Vector<double> > mem;
+ SolverCG<Vector<double> > solver (control, mem);
+
+ //SparseDirectUMFPACK direct_solver;
+ //direct_solver.initialize(filtered_A);
+
+ solver.solve(filtered_A,sol_total,rhs,precon);
+
+ vspace.get_hang().distribute(sol_total);
+
+ //sol_total += sol_increment;
+
+}
+
--- /dev/null
+#ifndef LOCAL_ASSEMBLE_STANDARD
+#define LOCAL_ASSEMBLE_STANDARD
+
+#include "local_assemble_base.h"
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+
+#include <fstream>
+#include <iostream>
+#include <base/parameter_handler.h>
+
+
+template <int dim, typename DH=MGDoFHandler<dim> >
+ class LocalAssemble : public LocalAssembleBase<dim, DH>
+{
+public:
+#if deal_II_dimension != 1
+ LocalAssemble();
+
+ /** Given two face finite elements, assemble the two matrices.*/
+ virtual void assemble_face_terms (FEFaceValuesBase<dim> &fe_v,
+ FEFaceValuesBase<dim> &fe_n_v,
+ FullMatrix<double> &,
+ FullMatrix<double> &);
+
+ virtual void assemble_face_term
+ (const typename DH::active_cell_iterator& ,
+ const unsigned int,
+ const typename DH::active_cell_iterator& ,
+ const unsigned int,
+ FullMatrix<double> &,
+ FullMatrix<double> &);
+
+ virtual void assemble_face_term
+ (const typename DH::active_cell_iterator& ,
+ const unsigned int,
+ const unsigned int,
+ const typename DH::active_cell_iterator& ,
+ const unsigned int,
+ FullMatrix<double> &,
+ FullMatrix<double> &);
+
+ virtual void assemble_face_term
+ (const typename DH::active_cell_iterator& ,
+ const unsigned int,
+ const typename DH::active_cell_iterator& ,
+ const unsigned int,
+ const unsigned int,
+ FullMatrix<double> &,
+ FullMatrix<double> &);
+
+ /** A pointer to fe_face_values objects. */
+ SmartPointer<FEFaceValues<dim> > fe_face_v;
+
+ /** A pointer to fe_face_values objects. */
+ SmartPointer<FEFaceValues<dim> > fe_face_n_v;
+
+ /** A pointer to fe_face_values objects. */
+ SmartPointer<FESubfaceValues<dim> > fe_sub_face_v;
+#endif
+
+};
+
+#endif
--- /dev/null
+#include "local_assemble.h"
+#include <ostream>
+
+#if deal_II_dimension != 1
+
+template <int dim, typename DH>
+LocalAssemble<dim,DH>::LocalAssemble() :
+ fe_face_v(0, "Local Assemble FeFaceValues Pointer"),
+ fe_face_n_v(0, "Local Assemble FeFaceValues on Neighbor Pointer"),
+ fe_sub_face_v(0, "Local Assemble FesubfaceValues Pointer")
+{
+}
+
+template <int dim, typename DH>
+void LocalAssemble<dim,DH>::assemble_face_terms (FEFaceValuesBase<dim> &,
+ FEFaceValuesBase<dim> &,
+ FullMatrix<double> &,
+ FullMatrix<double> &)
+{
+ Assert(false, ExcPureFunctionCalled());
+}
+
+
+template <int dim, typename DH>
+void LocalAssemble<dim,DH>::assemble_face_term
+(const typename DH::active_cell_iterator& cell,
+ const unsigned int face_no,
+ const typename DH::active_cell_iterator& n_cell,
+ const unsigned int n_face_no,
+ FullMatrix<double> &cell_m,
+ FullMatrix<double> &cell_n_m)
+{
+ this->fe_face_v->reinit(cell, face_no);
+ this->fe_face_n_v->reinit(n_cell, n_face_no);
+ assemble_face_terms(*fe_face_v, *fe_face_n_v, cell_m, cell_n_m);
+}
+
+template <int dim, typename DH>
+void LocalAssemble<dim,DH>::assemble_face_term
+(const typename DH::active_cell_iterator& cell,
+ const unsigned int face_no,
+ const typename DH::active_cell_iterator& n_cell,
+ const unsigned int n_face_no,
+ const unsigned int n_sub_face_no,
+ FullMatrix<double> &cell_m,
+ FullMatrix<double> &cell_n_m)
+{
+ fe_face_v->reinit(cell, face_no);
+ fe_sub_face_v->reinit(n_cell, n_face_no, n_sub_face_no);
+ assemble_face_terms(*fe_face_v, *fe_sub_face_v, cell_m, cell_n_m);
+}
+
+template <int dim, typename DH>
+void LocalAssemble<dim,DH>::assemble_face_term
+(const typename DH::active_cell_iterator& cell,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const typename DH::active_cell_iterator& n_cell,
+ const unsigned int n_face_no,
+ FullMatrix<double> &cell_m,
+ FullMatrix<double> &cell_n_m)
+{
+ fe_sub_face_v->reinit(cell, face_no, subface_no);
+ fe_face_v->reinit(n_cell, n_face_no);
+ assemble_face_terms(*fe_sub_face_v, *fe_face_v, cell_m, cell_n_m);
+}
+#endif
--- /dev/null
+#ifndef __deal2__local_assemble_base_h
+#define __deal2__local_assemble_base_h
+
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <fe/fe_values.h>
+#include "assemble_flags.h"
+
+#include <fstream>
+#include <iostream>
+
+using namespace dealii;
+
+/** Base class for local assemblers. This is the object you should
+ derive your own classes from, if you want to use the VirtualMatrix
+ class or the MyTools::assemble utility. It provides a common
+ interface to the procedure for assembling locally the matrices. All
+ functions are virtual and would throw an exception if called. You
+ should really implement your own class deriving it from this one. An
+ example is provided in the LocalAssembleLaplace class, which
+ assemble the Laplace equations locally for continuous finite
+ elements.
+
+ All methods work basically in the same way. They take a FEValues
+ object and a FullMatrix one as input, and fill the FullMatrix with
+ appropriate values. Exceptions are the face and boundary terms,
+ which take FEFaceValues as input.
+
+ This object is supposed to be used in a Virtual Matrix via the
+ VirtualMatrix::enter method.
+
+ Here the only thing which is actually done is to create the internal
+ AssembleFlags object upon construction, which is used to determine which
+ of the methods will be called by the virtual matrix.
+ */
+template <int dim, typename DH=MGDoFHandler<dim> >
+class LocalAssembleBase : public Subscriptor
+{
+ public:
+ virtual ~LocalAssembleBase() {};
+ /** This object will be called for each cell of a triangulation. */
+ // virtual void assemble_cell_term(const FEValues<dim>& fe_v,
+ // FullMatrix<double> &u_v_matrix) const;
+ virtual void assemble_cell_term
+ (const typename DH::active_cell_iterator&,
+ FullMatrix<double> &cell_m);
+
+ /** This object will be called for each boundary face of a triangulation.*/
+ virtual void assemble_boundary_term
+ (const typename DH::active_cell_iterator&, const unsigned int,
+ FullMatrix<double> &cell_m);
+
+ /** This object will be called for each face of a
+ triangulation. This one is called when the face is shared
+ by neighbors on the same level of refinement. */
+ virtual void assemble_face_term
+ (const typename DH::active_cell_iterator &cell,
+ const unsigned int face_no,
+ const typename DH::active_cell_iterator &neighbor,
+ const unsigned int neighbor_face_no,
+ FullMatrix<double> &cell_m,
+ FullMatrix<double> &neighbor_cell_m);
+
+ /** This object will be called for each face of a
+ triangulation. This one is called when the current face is
+ coarser. */
+ virtual void assemble_face_term
+ (const typename DH::active_cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int sub_face_no,
+ const typename DH::active_cell_iterator &neighbor,
+ const unsigned int neighbor_face_no,
+ FullMatrix<double> &cell_m,
+ FullMatrix<double> &neighbor_cell_m);
+
+ /** This object will be called for each face of a
+ triangulation. This one is called when the current face is
+ finer. */
+ virtual void assemble_face_term
+ (const typename DH::active_cell_iterator &cell,
+ const unsigned int face_no,
+ const typename DH::active_cell_iterator &neighbor,
+ const unsigned int neighbor_face_no,
+ const unsigned int neighbor_subface_no,
+ FullMatrix<double> &cell_m,
+ FullMatrix<double> &neighbor_cell_m);
+
+ /** Assemble rhs. This will be called for each cell.*/
+ virtual void assemble_rhs_term
+ (const typename DH::active_cell_iterator&,
+ Vector<double> &);
+
+ /** This object will be called for each boundary face of a triangulation.*/
+ virtual void assemble_rhs_boundary_term
+ (const typename DH::active_cell_iterator&, const unsigned int,
+ Vector<double> &);
+
+ AssembleFlags flags;
+};
+
+#endif
--- /dev/null
+#include "local_assemble_base.h"
+
+template <int dim, typename DH>
+void LocalAssembleBase<dim,DH>::assemble_cell_term (const typename DH::active_cell_iterator&,
+ FullMatrix<double> &)
+{
+ AssertThrow(false, ExcPureFunctionCalled());
+}
+
+template <int dim, typename DH>
+void LocalAssembleBase<dim,DH>::assemble_boundary_term
+(const typename DH::active_cell_iterator&,const unsigned int,
+ FullMatrix<double> &)
+{
+ AssertThrow(false, ExcPureFunctionCalled());
+}
+
+template <int dim, typename DH>
+void LocalAssembleBase<dim,DH>::assemble_face_term
+(const typename DH::active_cell_iterator&,const unsigned int,
+ const typename DH::active_cell_iterator&,const unsigned int,
+ FullMatrix<double> &,
+ FullMatrix<double> &)
+{
+ AssertThrow(false, ExcPureFunctionCalled());
+}
+
+template <int dim, typename DH>
+void LocalAssembleBase<dim,DH>::assemble_face_term
+(const typename DH::active_cell_iterator&,
+ const unsigned int, const unsigned int,
+ const typename DH::active_cell_iterator&,const unsigned int,
+ FullMatrix<double> &,
+ FullMatrix<double> &)
+{
+ AssertThrow(false, ExcPureFunctionCalled());
+}
+
+template <int dim, typename DH>
+void LocalAssembleBase<dim,DH>::assemble_face_term
+(const typename DH::active_cell_iterator&,const unsigned int,
+ const typename DH::active_cell_iterator&,
+ const unsigned int,
+ const unsigned int,
+ FullMatrix<double> &,
+ FullMatrix<double> &)
+{
+ AssertThrow(false, ExcPureFunctionCalled());
+}
+
+
+template <int dim, typename DH>
+void LocalAssembleBase<dim,DH>::assemble_rhs_term
+(const typename DH::active_cell_iterator&,
+ Vector<double> &)
+{
+ AssertThrow(false, ExcPureFunctionCalled());
+}
+
+template <int dim, typename DH>
+void LocalAssembleBase<dim,DH>::assemble_rhs_boundary_term
+(const typename DH::active_cell_iterator&,
+ const unsigned,
+ Vector<double> &)
+{
+ AssertThrow(false, ExcPureFunctionCalled());
+}
+
--- /dev/null
+#ifndef LOCAL_ASSEMBLE_ELASTIC_MATRIX
+#define LOCAL_ASSEMBLE_ELASTIC_MATRIX
+
+#include "local_assemble.h"
+#include "parsed_symmetric_tensor_function.h"
+#include <base/parsed_function.h>
+
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+
+#include <fstream>
+#include <iostream>
+#include <base/parameter_handler.h>
+
+
+template <int dim>
+class LocalAssembleElasticMatrix : public LocalAssembleBase<dim>
+{
+ public:
+
+ LocalAssembleElasticMatrix();
+
+ ~LocalAssembleElasticMatrix();
+
+ virtual void assemble_cell_term
+ (const typename MGDoFHandler<dim>::active_cell_iterator& cell,
+ FullMatrix<double> &);
+
+
+ void reinit(FiniteElement<dim>&);
+
+ void parameters(ParameterHandler &prm);
+
+ private:
+ /** A pointer to the finite element.*/
+ SmartPointer<FiniteElement<dim> > fe;
+
+ /** A pointer to the rhs function.*/
+ ParsedSymmetricTensorFunction<4, dim> C;
+
+ /** A pointer to fe_values objects. */
+ SmartPointer<FEValues<dim> > fe_v;
+
+};
+
+#endif
--- /dev/null
+#include "../include/local_assemble_elastic_matrix.h"
+#include <base/quadrature_lib.h>
+
+template <int dim>
+LocalAssembleElasticMatrix<dim>::LocalAssembleElasticMatrix() :
+fe(0, "Local Assemble Fe Pointer"),
+ fe_v(0, "Local Assemble FeValues Pointer")
+{
+}
+
+
+template <int dim>
+void LocalAssembleElasticMatrix<dim>::reinit (FiniteElement<dim>& myfe)
+{
+ if(fe_v) {
+ FEValues<dim> * p = fe_v;
+ fe_v = 0;
+ delete p;
+ }
+ fe = &myfe;
+ QGauss<dim> quadrature(2*fe->degree + 1);
+ QGauss<dim-1> face_quadrature(2*fe->degree + 1);
+ UpdateFlags flags (update_values |
+ update_gradients |
+ update_q_points |
+ update_JxW_values);
+ fe_v = new FEValues<dim>(*fe, quadrature, flags);
+
+ this->flags = assemble_cell;
+
+}
+
+template <int dim>
+LocalAssembleElasticMatrix<dim>::~LocalAssembleElasticMatrix()
+{
+ if(fe_v) {
+ FEValues<dim> * p = fe_v;
+ fe_v = 0;
+ delete p;
+ }
+ fe = 0;
+}
+
+template <int dim>
+void LocalAssembleElasticMatrix<dim>::parameters(ParameterHandler &prm)
+{
+ prm.enter_subsection("Elastic Moduli");
+ C.parse_parameters(prm);
+ prm.leave_subsection();
+}
+
+template <int dim>
+void LocalAssembleElasticMatrix<dim>::assemble_cell_term
+(const typename MGDoFHandler<dim>::active_cell_iterator& cell,
+ FullMatrix<double> &cell_m)
+{
+ cell_m = 0;
+ Assert(fe, ExcNotInitialized());
+
+ fe_v->reinit(cell);
+
+ std::vector<Point<dim> > points(fe_v->n_quadrature_points);
+ points = fe_v->get_quadrature_points();
+ unsigned int comp_i, comp_j;
+
+ for(unsigned int point =0; point<fe_v->n_quadrature_points; ++point) {
+ SymmetricTensor<4,dim> C_qp = C(points[point]);
+
+ for (unsigned int i=0; i<fe_v->dofs_per_cell; ++i) {
+ comp_i = fe->system_to_component_index(i).first;
+ for (unsigned int j=0; j<fe_v->dofs_per_cell; ++j) {
+ comp_j = fe->system_to_component_index(j).first;
+
+ for(unsigned int b=0; b<dim; ++b) {
+ for(unsigned int n=0; n<dim; ++n) {
+ cell_m(i,j) += C_qp[comp_i][b][comp_j][n] *
+ fe_v->shape_grad(j, point)[n] *
+ fe_v->shape_grad(i, point)[b] *
+ fe_v->JxW(point);
+ } //n
+ }//b
+ } //j
+ }//i
+ }//qp
+}
--- /dev/null
+#ifndef LOCAL_ASSEMBLE_ELASTIC_RHS
+#define LOCAL_ASSEMBLE_ELASTIC_RHS
+
+#include "local_assemble.h"
+#include "parsed_symmetric_tensor_function.h"
+#include <base/parsed_function.h>
+
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+
+#include <fstream>
+#include <iostream>
+#include <base/parameter_handler.h>
+
+
+template <int dim>
+class LocalAssembleElasticRHS : public LocalAssembleBase<dim>
+{
+ public:
+
+ LocalAssembleElasticRHS();
+
+ ~LocalAssembleElasticRHS();
+
+ virtual void assemble_rhs_term
+ (const typename MGDoFHandler<dim>::active_cell_iterator&,
+ Vector<double> &);
+
+ /** This object will be called for each boundary face of a triangulation.*/
+ virtual void assemble_rhs_boundary_term
+ (const typename MGDoFHandler<dim>::active_cell_iterator&, const unsigned int,
+ Vector<double> &);
+
+
+ void reinit(FiniteElement<dim>&,
+ Function<dim> &,
+ Function<dim> &,
+ std::map<char, std::vector<bool> > &);
+
+ private:
+ /** A pointer to the finite element.*/
+ SmartPointer<FiniteElement<dim> > fe;
+
+ /** A pointer to the rhs function.*/
+ SmartPointer<Function<dim> > rhs;
+
+ /** A pointer to the Neumann function.*/
+ SmartPointer<Function<dim> > neumann;
+
+ /** A map of ids and components for neumann bc. */
+ std::map<char, std::vector<bool> > neumann_map;
+
+ /** A pointer to fe_values objects. */
+ SmartPointer<FEValues<dim> > fe_v;
+
+ /** A pointer to fe_face_values objects. */
+ SmartPointer<FEFaceValues<dim> > fe_face_v;
+
+};
+
+#endif
--- /dev/null
+#include "../include/local_assemble_elastic_rhs.h"
+#include <base/quadrature_lib.h>
+
+template <int dim>
+LocalAssembleElasticRHS<dim>::LocalAssembleElasticRHS() :
+fe(0, "Local Assemble Fe Pointer"),
+ fe_v(0, "Local Assemble FeValues Pointer"),
+ fe_face_v(0, "Local Assemble FeFaceValues Pointer")
+{
+}
+
+
+template <int dim>
+void LocalAssembleElasticRHS<dim>::reinit (FiniteElement<dim> &myfe,
+ Function<dim> &bf,
+ Function<dim> &nbc,
+ std::map<char, std::vector<bool> > & n_map)
+{
+ if(fe_v) {
+ FEValues<dim> * p = fe_v;
+ fe_v = 0;
+ delete p;
+ }
+ if(fe_face_v) {
+ FEFaceValues<dim> * p = fe_face_v;
+ fe_face_v = 0;
+ delete p;
+ }
+ fe = &myfe;
+ QGauss<dim> quadrature(2*fe->degree + 1);
+ QGauss<dim-1> face_quadrature(2*fe->degree + 1);
+ UpdateFlags flags (update_values |
+ update_gradients |
+ update_q_points |
+ update_JxW_values);
+ fe_v = new FEValues<dim>(*fe, quadrature, flags);
+
+ fe_face_v = new FEFaceValues<dim>(*fe, face_quadrature, flags);
+
+ this->flags =
+ assemble_rhs_cell|
+ assemble_rhs_boundary;
+
+ rhs = &bf;
+ neumann = &nbc;
+ neumann_map = n_map;
+}
+
+template <int dim>
+LocalAssembleElasticRHS<dim>::~LocalAssembleElasticRHS()
+{
+ if(fe_v) {
+ FEValues<dim> * p = fe_v;
+ fe_v = 0;
+ delete p;
+ }
+ if(fe_face_v) {
+ FEFaceValues<dim> * p = fe_face_v;
+ fe_face_v = 0;
+ delete p;
+ }
+ fe = 0;
+}
+
+
+template <int dim>
+void LocalAssembleElasticRHS<dim>::assemble_rhs_term
+(const typename MGDoFHandler<dim>::active_cell_iterator& cell,
+ Vector<double> &cell_rhs)
+{
+ cell_rhs = 0;
+ Assert(fe, ExcNotInitialized());
+
+ fe_v->reinit(cell);
+
+ unsigned int size = fe->n_components();
+ unsigned int n_q_points = fe_v->n_quadrature_points;
+
+ std::vector<Vector<double> > load_vector (n_q_points, Vector<double>(size) );
+
+ /* Evaluate rhs and solution on the quadrature points. */
+ rhs->vector_value_list (fe_v->get_quadrature_points(), load_vector);
+
+ unsigned int comp_i;
+ for (unsigned int i=0; i<fe_v->dofs_per_cell; ++i) {
+ comp_i = fe->system_to_component_index(i).first;
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point) {
+ cell_rhs(i) +=
+ ( load_vector[q_point](comp_i) *
+ fe_v->shape_value(i,q_point) *
+ fe_v->JxW(q_point) );
+ }
+ }
+}
+
+template<int dim>
+void LocalAssembleElasticRHS<dim>::assemble_rhs_boundary_term
+(const typename MGDoFHandler<dim>::active_cell_iterator& cell,
+ const unsigned int face_no,
+ Vector<double> &cell_rhs)
+{
+ cell_rhs = 0;
+ // See if we need to do anything here
+ char id = cell->face(face_no)->boundary_indicator();
+ if(neumann_map.find(id) == neumann_map.end())
+ return;
+
+ std::vector<bool> & filter = neumann_map[id];
+
+ Assert(fe, ExcNotInitialized());
+ fe_face_v->reinit(cell,face_no);
+
+ unsigned int size = fe->n_components();
+ unsigned int n_q_points = fe_face_v->n_quadrature_points;
+
+ /** Vector of boundary values.*/
+ std::vector<Vector<double> > neumann_vector (n_q_points, Vector<double>(size) );
+
+ /* Evaluate rhs and solution on the quadrature points. */
+ neumann->vector_value_list (fe_face_v->get_quadrature_points(), neumann_vector);
+
+ unsigned int comp_i;
+ for (unsigned int i=0; i<fe_face_v->dofs_per_cell; ++i) {
+ comp_i = fe->system_to_component_index(i).first;
+ if(filter[comp_i])
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point) {
+ cell_rhs(i) +=
+ ( neumann_vector[q_point](comp_i) *
+ fe_face_v->shape_value(i,q_point) *
+ fe_face_v->JxW(q_point) );
+ }
+ }
+}
--- /dev/null
+#ifndef LOCAL_ASSEMBLE_HARDENING_PROJECT
+#define LOCAL_ASSEMBLE_HARDENING_PROJECT
+
+#include "local_assemble.h"
+#include "parsed_symmetric_tensor_function.h"
+#include <base/parsed_function.h>
+
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+
+#include <fstream>
+#include <iostream>
+#include <base/parameter_handler.h>
+
+
+template <int dim>
+class LocalAssemblePlasticProject : public LocalAssembleBase<dim>
+{
+ public:
+
+ LocalAssemblePlasticProject();
+
+ ~LocalAssemblePlasticProject();
+
+ virtual void assemble_rhs_term
+ (const typename MGDoFHandler<dim>::active_cell_iterator&,
+ Vector<double> &);
+
+ void parameters(ParameterHandler &prm);
+
+
+ void reinit(FiniteElement<dim>&,
+ Table<2, SymmetricTensor<2,dim> > &);
+
+ private:
+ /** A pointer to the finite element.*/
+ SmartPointer<FiniteElement<dim> > fe;
+
+ /** A pointer to the plastic strain */
+ Table<dim,SymmetricTensor<2,dim> > plastic_strain;
+
+ /** A pointer to fe_values objects. */
+ SmartPointer<FEValues<dim> > fe_v;
+
+ ParsedSymmetricTensorFunction<4, dim> C;
+
+};
+
+#endif
--- /dev/null
+#ifndef LOCAL_ASSEMBLE_MASS_MATRIX
+#define LOCAL_ASSEMBLE_MASS_MATRIX
+
+#include "local_assemble.h"
+
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+
+template <int dim, typename DH>
+class LocalAssembleMass :
+public LocalAssembleBase<dim, DH>
+{
+ public:
+
+ LocalAssembleMass();
+
+ ~LocalAssembleMass();
+
+ virtual void assemble_cell_term
+ (const typename DH::active_cell_iterator& cell,
+ FullMatrix<double> &);
+
+ virtual void assemble_rhs_term
+ (const typename DH::active_cell_iterator&,
+ Vector<double> &);
+
+ void reinit(FiniteElement<dim>&,
+ Function<dim> &);
+
+ private:
+ /** A pointer to the finite element.*/
+ SmartPointer<FiniteElement<dim> > fe;
+
+ /** A pointer to the rhs function.*/
+ SmartPointer<Function<dim> > rhs;
+
+ /** A pointer to fe_values objects. */
+ SmartPointer<FEValues<dim> > fe_v;
+};
+
+#endif
--- /dev/null
+#include "../include/local_assemble_mass.h"
+#include "utilities.h"
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+
+template <int dim, typename DH>
+LocalAssembleMass<dim, DH>::LocalAssembleMass() :
+ fe(0, "Local Assemble Fe Pointer"),
+ fe_v(0, "Local Assemble FeValues Pointer")
+{
+}
+
+
+template <int dim, typename DH>
+void LocalAssembleMass<dim, DH>::reinit (FiniteElement<dim> &myfe,
+ Function<dim> &f)
+{
+ smart_delete(fe_v);
+ fe = &myfe;
+ QGauss<dim> quadrature(2*fe->degree + 1);
+ UpdateFlags flags (update_values |
+ update_q_points |
+ update_JxW_values);
+ fe_v = new FEValues<dim>(*fe, quadrature, flags);
+
+ this->flags =
+ assemble_cell|
+ assemble_rhs_cell;
+
+ rhs = &f;
+}
+
+template <int dim, typename DH>
+LocalAssembleMass<dim, DH>::~LocalAssembleMass()
+{
+ smart_delete(fe_v);
+ fe = 0;
+}
+
+template <int dim, typename DH>
+void LocalAssembleMass<dim, DH>::assemble_cell_term
+(const typename DH::active_cell_iterator& cell,
+ FullMatrix<double> &cell_m)
+{
+ cell_m = 0;
+ Assert(fe, ExcNotInitialized());
+
+ fe_v->reinit(cell);
+ unsigned int comp_i = 0, comp_j = 0;
+ for (unsigned int i=0; i<fe_v->dofs_per_cell; ++i) {
+ comp_i = fe->system_to_component_index(i).first;
+ for (unsigned int j=0; j<fe_v->dofs_per_cell; ++j) {
+ comp_j = fe->system_to_component_index(j).first;
+ if(comp_i == comp_j)
+ for(unsigned int q_point =0; q_point<fe_v->n_quadrature_points; ++q_point) {
+ cell_m(i,j) += ( fe_v->shape_value(j,q_point) *
+ fe_v->shape_value(i,q_point) *
+ fe_v->JxW(q_point) );
+ }
+ }
+ }
+}
+
+
+template <int dim, typename DH>
+void LocalAssembleMass<dim, DH>::assemble_rhs_term
+(const typename DH::active_cell_iterator& cell,
+ Vector<double> &cell_rhs)
+{
+ cell_rhs = 0;
+ Assert(fe, ExcNotInitialized());
+
+ fe_v->reinit(cell);
+
+ unsigned int size = fe->n_components();
+ unsigned int n_q_points = fe_v->n_quadrature_points;
+
+ std::vector<Vector<double> > load_vector (n_q_points, Vector<double>(size) );
+
+ /* Evaluate rhs and solution on the quadrature points. */
+ rhs->vector_value_list (fe_v->get_quadrature_points(), load_vector);
+
+ unsigned int comp_i;
+ for (unsigned int i=0; i<fe_v->dofs_per_cell; ++i) {
+ comp_i = fe->system_to_component_index(i).first;
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point) {
+ cell_rhs(i) +=
+ ( load_vector[q_point](comp_i) *
+ fe_v->shape_value(i,q_point) *
+ fe_v->JxW(q_point) );
+ }
+ }
+}
--- /dev/null
+#ifndef LOCAL_ASSEMBLE_PLASTIC_PROJECT
+#define LOCAL_ASSEMBLE_PLASTIC_PROJECT
+
+#include "local_assemble.h"
+#include "parsed_symmetric_tensor_function.h"
+#include <base/parsed_function.h>
+
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+
+#include <fstream>
+#include <iostream>
+#include <base/parameter_handler.h>
+
+
+template <int dim>
+class LocalAssemblePlasticProject : public LocalAssembleBase<dim>
+{
+ public:
+
+ LocalAssemblePlasticProject();
+
+ ~LocalAssemblePlasticProject();
+
+ virtual void assemble_rhs_term
+ (const typename MGDoFHandler<dim>::active_cell_iterator&,
+ Vector<double> &);
+
+ void parameters(ParameterHandler &prm);
+
+
+ void reinit(FiniteElement<dim>&,
+ Table<2, SymmetricTensor<2,dim> > &);
+
+ private:
+ /** A pointer to the finite element.*/
+ SmartPointer<FiniteElement<dim> > fe;
+
+ /** A pointer to the plastic strain */
+ Table<dim,SymmetricTensor<2,dim> > plastic_strain;
+
+ /** A pointer to fe_values objects. */
+ SmartPointer<FEValues<dim> > fe_v;
+
+ ParsedSymmetricTensorFunction<4, dim> C;
+
+};
+
+#endif
--- /dev/null
+#include "../include/local_assemble_plastic_project.h"
+#include <base/quadrature_lib.h>
+
+template <int dim>
+LocalAssemblePlasticProject<dim>::LocalAssemblePlasticProject() :
+fe(0, "Local Assemble Fe Pointer"),
+ fe_v(0, "Local Assemble FeValues Pointer")
+{
+}
+
+
+template <int dim>
+void LocalAssemblePlasticProject<dim>::reinit (FiniteElement<dim> &myfe,
+ Table<2, SymmetricTensor<2,dim> > &my_plastic_strain)
+{
+ if(fe_v) {
+ FEValues<dim> * p = fe_v;
+ fe_v = 0;
+ delete p;
+ }
+ fe = &myfe;
+ QGauss<dim> quadrature(2*fe->degree + 1);
+ QGauss<dim-1> face_quadrature(2*fe->degree + 1);
+ UpdateFlags flags (update_gradients |
+ update_q_points |
+ update_JxW_values);
+ fe_v = new FEValues<dim>(*fe, quadrature, flags);
+
+ this->flags = assemble_rhs_cell;
+
+ plastic_strain = my_plastic_strain;
+}
+
+template <int dim>
+LocalAssemblePlasticProject<dim>::~LocalAssemblePlasticProject()
+{
+ if(fe_v) {
+ FEValues<dim> * p = fe_v;
+ fe_v = 0;
+ delete p;
+ }
+ fe = 0;
+}
+
+template <int dim>
+void LocalAssemblePlasticProject<dim>::parameters(ParameterHandler &prm)
+{
+ prm.enter_subsection("Elastic Moduli");
+ C.parse_parameters(prm);
+ prm.leave_subsection();
+}
+
+
+template <int dim>
+void LocalAssemblePlasticProject<dim>::assemble_rhs_term
+(const typename MGDoFHandler<dim>::active_cell_iterator& cell,
+ Vector<double> &cell_rhs)
+{
+ cell_rhs = 0;
+ Assert(fe, ExcNotInitialized());
+
+ fe_v->reinit(cell);
+
+ unsigned int size = fe->n_components();
+ unsigned int n_q_points = fe_v->n_quadrature_points;
+
+ std::vector<Vector<double> > load_vector (n_q_points, Vector<double>(size) );
+ std::vector<Point<dim> > points = fe_v->get_quadrature_points();
+
+ for (unsigned int qp=0; qp<n_q_points; ++qp) {
+
+ SymmetricTensor<4,dim> C_qp = C(points[qp]);
+
+ for (unsigned int i=0; i<fe_v->dofs_per_cell; ++i) {
+ unsigned int comp_i = fe->system_to_component_index(i).first;
+
+ for(unsigned int b=0; b<dim; ++b) {
+ for(unsigned int m=0; m<dim; ++m) {
+ for(unsigned int n=0; n<dim; ++n) {
+
+ cell_rhs(i) += ( C_qp[comp_i][b][m][n] *
+ plastic_strain(cell->index(), qp)[m][n] *
+ fe_v->shape_grad(i, qp)[b] *
+ fe_v->JxW(qp) );
+
+ } //n
+ } //m
+ } //b
+ } //i
+ } //qp
+}
+
+
--- /dev/null
+#ifndef LOCAL_ASSEMBLE_SCALAR_PROJECT
+#define LOCAL_ASSEMBLE_SCALAR_PROJECT
+
+#include "local_assemble.h"
+#include <base/parsed_function.h>
+
+#include <base/logstream.h>
+#include <base/smartpointer.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+
+#include <fstream>
+#include <iostream>
+#include <base/parameter_handler.h>
+
+
+template <int dim>
+class LocalAssembleScalarProject : public LocalAssembleBase<dim>
+{
+ public:
+
+ LocalAssembleScalarProject();
+
+ ~LocalAssembleScalarProject();
+
+ virtual void assemble_rhs_term
+ (const typename MGDoFHandler<dim>::active_cell_iterator&,
+ Vector<double> &);
+
+
+ void reinit(FiniteElement<dim>&,
+ Table<2, double > &,
+ Table<2, double > &);
+
+ private:
+ /** A pointer to the finite element.*/
+ SmartPointer<FiniteElement<dim> > fe;
+
+ /** A pointer to the plastic strain */
+ Table<dim,double > hard_table;
+ Table<dim,double > iter_table;
+
+ /** A pointer to fe_values objects. */
+ SmartPointer<FEValues<dim> > fe_v;
+
+
+};
+
+#endif
--- /dev/null
+#include "../include/local_assemble_scalar_project.h"
+#include <base/quadrature_lib.h>
+
+template <int dim>
+LocalAssembleScalarProject<dim>::LocalAssembleScalarProject() :
+fe(0, "Local Assemble Fe Pointer"),
+ fe_v(0, "Local Assemble FeValues Pointer")
+{
+}
+
+
+template <int dim>
+void LocalAssembleScalarProject<dim>::reinit (FiniteElement<dim> &myfe,
+ Table<2, double> &my_hard,
+ Table<2, double> &my_iter)
+{
+ if(fe_v) {
+ FEValues<dim> * p = fe_v;
+ fe_v = 0;
+ delete p;
+ }
+ fe = &myfe;
+ QGauss<dim> quadrature(2*fe->degree + 1);
+ QGauss<dim-1> face_quadrature(2*fe->degree + 1);
+ UpdateFlags flags (update_values |
+ update_q_points |
+ update_JxW_values);
+ fe_v = new FEValues<dim>(*fe, quadrature, flags);
+
+ this->flags = assemble_rhs_cell;
+
+ hard_table = my_hard;
+ iter_table = my_iter;
+}
+
+template <int dim>
+LocalAssembleScalarProject<dim>::~LocalAssembleScalarProject()
+{
+ if(fe_v) {
+ FEValues<dim> * p = fe_v;
+ fe_v = 0;
+ delete p;
+ }
+ fe = 0;
+}
+
+
+template <int dim>
+void LocalAssembleScalarProject<dim>::assemble_rhs_term
+(const typename MGDoFHandler<dim>::active_cell_iterator& cell,
+ Vector<double> &cell_rhs)
+{
+ cell_rhs = 0;
+ Assert(fe, ExcNotInitialized());
+
+ fe_v->reinit(cell);
+
+ unsigned int n_q_points = fe_v->n_quadrature_points;
+
+ for (unsigned int qp=0; qp<n_q_points; ++qp) {
+
+ for (unsigned int i=0; i<fe_v->dofs_per_cell; ++i) {
+ unsigned int comp_i = fe->system_to_component_index(i).first;
+
+ if(comp_i == 0) {
+
+ cell_rhs(i) += ( hard_table(cell->index(), qp) *
+ fe_v->shape_value(i, qp) *
+ fe_v->JxW(qp) );
+ }
+
+ if(comp_i == 1) {
+
+ cell_rhs(i) += ( iter_table(cell->index(), qp) *
+ fe_v->shape_value(i, qp) *
+ fe_v->JxW(qp) );
+ }
+
+
+ } //i
+ } //qp
+}
+
+
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id: fe.h,v 1.124 2005/09/17 09:19:19 guido Exp $
+// Version: $Name: $
+//
+// Copyright (C) 2005 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__my_tools_h
+#define __deal2__my_tools_h
+
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <multigrid/mg_base.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <fe/fe.h>
+#include <base/quadrature_lib.h>
+
+#include <dofs/dof_constraints.h>
+#include "assemble_flags.h"
+#include "local_assemble.h"
+#include "vector_space.h"
+
+/**
+ * Some assembly routines. These are an extension of the already
+ * existing deal.II assembly routines that make use of the support
+ * class LocalAssembleBase. These routines call functions of the local
+ * assemblers passing local matrices or vectors and the cell. It is
+ * responsability of the local assemblers to fill in the informations
+ * correctly. Unlike step-12 of the deal.II library, the cell iterator
+ * is passed, not the fe_values. This means that the local assembler
+ * should build its own fe_values object. It is its responsability to
+ * select all the outer details on how to actually compute these
+ * matrices...
+ *
+ *
+ * @author Luca Heltai,
+ * 2005, 2008
+ */
+
+namespace MyTools {
+ using namespace std;
+
+ /** DG-Like assembly routine for matrices. Hanging node
+ * constraints are not taken into account, and one should do this
+ * outside. However face terms are correctly handled, and this
+ * method is capable of assemblying fully DG objects. */
+ template <int dim, typename DH, typename MATRIX>
+ void assemble(DH &dof_handler, MATRIX& system_matrix,
+ LocalAssembleBase<dim, DH> &local,
+ const unsigned int this_mpi_process = 0)
+ {
+ const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+ std::vector<unsigned int> dofs (dofs_per_cell);
+ std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
+
+ // Now we create the cell matrices and vectors. Here we need two
+ // cell matrices, both for face terms that include test functions
+ // <code>vi</code> (internal shape functions, i.e. shape functions
+ // of the current cell).
+ //
+ // i stands for internal, e stands for external
+ FullMatrix<double> ui_vi_matrix (dofs_per_cell, dofs_per_cell);
+ FullMatrix<double> ui_ve_matrix (dofs_per_cell, dofs_per_cell);
+
+ // Furthermore we need some cell iterators.
+ typename DH::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ // Now we start the loop over all
+ // active cells.
+ for (;cell!=endc; ++cell)
+ if( cell->subdomain_id() == this_mpi_process )
+ {
+ // On each cell we need to reset the <code>ui_vi_matrix</code>
+ // and <code>cell_vector</code> to zero, before assembling the
+ // cell terms.
+ ui_vi_matrix = 0;
+
+ // and call the function that assembles the cell terms. The
+ // first argument is the <code>FEValues</code> that was
+ // previously reinit'ed on the current cell.
+ if( local.flags & assemble_cell )
+ local.assemble_cell_term(cell, ui_vi_matrix);
+
+ // As in previous examples the vector `dofs' includes the
+ // dof_indices.
+ cell->get_dof_indices (dofs);
+
+ // This is the start of the nested loop over all faces.
+ if( local.flags & (assemble_boundary|assemble_face) ) {
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ {
+ // First we set the face iterator
+ typename DH::face_iterator face=cell->face(face_no);
+
+ // and clear the <code>ui_ve_matrix</code> on each face.
+ ui_ve_matrix = 0;
+
+ // Now we distinguish the four different cases in the
+ // ordering mentioned above. We start with faces belonging
+ // to the boundary of the domain.
+ if (face->at_boundary())
+ {
+ // and assemble the corresponding face terms.
+ if(local.flags & assemble_boundary)
+ local.assemble_boundary_term(cell, face_no,
+ ui_vi_matrix);
+ }
+ else if(local.flags & assemble_face)
+ {
+ // Now we are not on the boundary of the domain,
+ // therefore there must exist a neighboring cell.
+ typename DH::cell_iterator neighbor=
+ cell->neighbor(face_no);
+
+ // We proceed with the second and most complicated case:
+ // the neighboring cell is more refined than the current
+ // cell. As in deal.II neighboring cells are restricted
+ // to have a level difference of not more than one, the
+ // neighboring cell is known to be at most ONCE more
+ // refined than the current cell. Furthermore also the
+ // face is more refined, i.e. it has children. Here we
+ // note that the following part of code will not work
+ // for <code>dim==1</code>.
+ if (face->has_children())
+ {
+ // First we store which number the current cell has
+ // in the list of neighbors of the neighboring
+ // cell. Hence,
+ // neighbor-@>neighbor(neighbor_face_no) equals the
+ // current cell <code>cell</code>.
+ const unsigned int neighbor_face_no=
+ cell->neighbor_of_neighbor(face_no);
+
+
+ // We loop over subfaces
+ for (unsigned int subface_no=0;
+ subface_no<face->n_children(); ++subface_no)
+ {
+ // and set the cell iterator
+ // <code>neighbor_child</code> to the cell
+ // placed `behind' the current subface.
+ typename DH::active_cell_iterator
+ neighbor_child
+ = cell->neighbor_child_on_subface (face_no, subface_no);
+
+ // As these are quite complicated indirections
+ // which one does not usually get right at first
+ // attempt we check for the internal
+ // consistency.
+ Assert (neighbor_child->face(neighbor_face_no) == face->child(subface_no),
+ ExcInternalError());
+ Assert (!neighbor_child->has_children(), ExcInternalError());
+
+ // We need to reset the
+ // <code>ui_ve_matrix</code> on each subface
+ // because on each subface the <code>un</code>
+ // belong to different neighboring cells.
+ ui_ve_matrix = 0;
+
+ // As already mentioned above for the current
+ // case (case 2) we employ the
+ // <code>FESubfaceValues</code> of the current
+ // cell (here reinited for the current cell,
+ // face and subface) and we employ the
+ // FEFaceValues of the neighboring child cell.
+ local.assemble_face_term(cell, face_no, subface_no,
+ neighbor_child, neighbor_face_no,
+ ui_vi_matrix,
+ ui_ve_matrix);
+
+ // Then we get the dof indices of the
+ // neighbor_child cell
+ neighbor_child->get_dof_indices (dofs_neighbor);
+
+ // and distribute <code>ui_ve_matrix</code> to
+ // the system_matrix
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ system_matrix.add(dofs_neighbor[i], dofs[k],
+ ui_ve_matrix(i,k));
+ }
+ // End of <code>if (face-@>has_children())</code>
+ }
+ else
+ {
+ // We proceed with case 3, i.e. neighboring cell is
+ // of the same refinement level as the current cell.
+ if (neighbor->level() == cell->level())
+ {
+ // Like before we store which number the current
+ // cell has in the list of neighbors of the
+ // neighboring cell.
+ const unsigned int neighbor_face_no=cell->neighbor_of_neighbor(face_no);
+
+ local.assemble_face_term(cell, face_no,
+ neighbor, neighbor_face_no,
+ ui_vi_matrix,
+ ui_ve_matrix);
+ // End of <code>if (neighbor-@>level() ==
+ // cell-@>level())</code>
+ }
+ else
+ {
+ // Finally we consider case 4. When the
+ // neighboring cell is not finer and not of the
+ // same refinement level as the current cell it
+ // must be coarser.
+ Assert(neighbor->level() < cell->level(), ExcInternalError());
+
+ // Find out the how many'th face_no and
+ // subface_no the current face is w.r.t. the
+ // neighboring cell.
+ const std::pair<unsigned int, unsigned int> faceno_subfaceno=
+ cell->neighbor_of_coarser_neighbor(face_no);
+ const unsigned int neighbor_face_no=faceno_subfaceno.first,
+ neighbor_subface_no=faceno_subfaceno.second;
+
+ Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
+ neighbor_subface_no)
+ == cell,
+ ExcInternalError());
+
+ local.assemble_face_term(cell, face_no,
+ neighbor,
+ neighbor_face_no,
+ neighbor_subface_no,
+ ui_vi_matrix,
+ ui_ve_matrix);
+ }
+
+ // Now we get the dof indices of the
+ // <code>neighbor</code> cell,
+ neighbor->get_dof_indices (dofs_neighbor);
+
+ // and distribute the
+ // <code>ui_ve_matrix</code>.
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ system_matrix.add(dofs_neighbor[i], dofs[k],
+ ui_ve_matrix(i,k));
+ }
+ // End of <code>face not at boundary</code>:
+ }
+ // End of loop over all faces:
+ }
+ }
+
+ // Finally we distribute the
+ // <code>ui_vi_matrix</code>
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i,j));
+ }
+ }
+
+
+ /** MG assembly routine for matrices. Hanging node constraints are
+ * taken into account. However face terms are NOT correctly
+ * handled, that is, this method is currently not capable of
+ * handling fluxex.
+
+ Different levels can be assembled. If the level is -1, then the
+ active level is used.
+ */
+ template <int dim, typename MATRIX>
+ void assemble(MGDoFHandler<dim> &mg_dof_handler, ConstraintMatrix &hang,
+ int level, MATRIX& m, LocalAssembleBase<dim, MGDoFHandler<dim> > &local,
+ const unsigned int this_mpi_process = 0)
+ {
+ AssertThrow( !(local.flags & assemble_face), ExcNotImplemented());
+ // Cell iterators.
+ typename MGDoFHandler<dim>::active_cell_iterator cell, endc;
+ if(level < 0) {
+ cell = mg_dof_handler.begin_active();
+ endc = mg_dof_handler.end();
+ } else {
+ cell = mg_dof_handler.begin_active(level);
+ endc = mg_dof_handler.end(level);
+ }
+
+ unsigned int dofs_per_cell =
+ mg_dof_handler.get_fe().dofs_per_cell;
+
+ // dofs and neighbor dofs.
+ vector<unsigned int> dofs(dofs_per_cell);
+ vector<unsigned int> neighbor_dofs(dofs_per_cell);
+
+ FullMatrix<double> cell_m(dofs_per_cell, dofs_per_cell);
+ FullMatrix<double> cell_neighbor_m(dofs_per_cell, dofs_per_cell);
+
+ for (; cell!=endc; ++cell)
+ if( cell->subdomain_id() == this_mpi_process )
+ {
+ // Reset local matrices.
+ cell_m = 0;
+ cell_neighbor_m = 0;
+
+ // Update indices
+ if(level<0)
+ cell->get_dof_indices(dofs);
+ else
+ cell->get_mg_dof_indices (dofs);
+
+ // Assemble cell term if needed
+ if( local.flags & assemble_cell ) {
+
+ //
+ local.assemble_cell_term(cell, cell_m);
+ }
+ for(unsigned int n=0; n<dofs.size(); ++n)
+ hang.distribute_local_to_global(cell_m, dofs, m);
+#if (DEAL_II_DIM > 1)
+ // Assemble boundary or face terms if needed
+ if( local.flags & (assemble_boundary|assemble_face) ) {
+
+ for (unsigned int face_no=0;
+ face_no<GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
+ {
+ typename MGDoFHandler<dim>::face_iterator face=cell->face(face_no);
+
+ if (face->at_boundary())
+ {
+ if(local.flags & assemble_boundary) {
+ cell_m = 0;
+ local.assemble_boundary_term(cell, face_no, cell_m);
+ hang.distribute_local_to_global(cell_m, dofs, m);
+ }
+ } // End of Boundary terms
+ else if(local.flags & assemble_face) {
+ typename MGDoFHandler<dim>::cell_iterator neighbor=
+ cell->neighbor(face_no);;
+
+ if (face->has_children()) // The neighbor is more refined
+ {
+ const unsigned int n_face_no=
+ cell->neighbor_of_neighbor(face_no);
+
+
+ for (unsigned int subface_no=0;
+ subface_no<face->n_children(); ++subface_no)
+ {
+ typename MGDoFHandler<dim>::active_cell_iterator
+ neighbor_child
+ = cell->neighbor_child_on_subface (face_no, subface_no);
+
+ Assert (neighbor_child->face(n_face_no) == face->child(subface_no),
+ ExcInternalError());
+ Assert (!neighbor_child->has_children(), ExcInternalError());
+
+ cell_m = 0;
+ cell_neighbor_m = 0;
+ local.assemble_face_term(cell, face_no, subface_no,
+ neighbor_child, n_face_no,
+ cell_m, cell_neighbor_m);
+ if(level<0)
+ neighbor_child->get_dof_indices (neighbor_dofs);
+ else
+ neighbor_child->get_mg_dof_indices (neighbor_dofs);
+
+ hang.distribute_local_to_global(cell_m, dofs, m);
+ // TBA Distribute neighbors!!!
+ // hang.distribute_local_to_global(cell_m, neighbor_dofs, dofs, m);
+ //
+
+
+ }
+ } // face has children
+ else
+ {
+ if (neighbor->level() == cell->level())
+ {
+ const unsigned int n_face_no=cell->neighbor_of_neighbor(face_no);
+
+ cell_neighbor_m = 0;
+ cell_m = 0;
+ local.assemble_face_term(cell, face_no,
+ neighbor, n_face_no,
+ cell_m,
+ cell_neighbor_m);
+
+ if(level<0)
+ neighbor->get_dof_indices (neighbor_dofs);
+ else
+ neighbor->get_mg_dof_indices (neighbor_dofs);
+
+ hang.distribute_local_to_global(cell_m, dofs, m);
+ // TBA: distribute neighbor dofs....
+ // hang.distribute_local_to_global(cell_m, neighbor_dofs, dofs, m);
+ //
+
+ } // Same level on the two neighbors
+ else
+ {
+ Assert(neighbor->level() < cell->level(), ExcInternalError());
+
+ const std::pair<unsigned int, unsigned int> faceno_subfaceno=
+ cell->neighbor_of_coarser_neighbor(face_no);
+ const unsigned int neighbor_face_no=faceno_subfaceno.first,
+ neighbor_subface_no=faceno_subfaceno.second;
+
+ Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
+ neighbor_subface_no)
+ == cell,
+ ExcInternalError());
+
+ cell_neighbor_m = 0;
+ cell_m = 0;
+ local.assemble_face_term(cell, face_no,
+ neighbor, neighbor_face_no,
+ neighbor_subface_no,
+ cell_m, cell_neighbor_m);
+
+ if(level<0)
+ neighbor->get_dof_indices (neighbor_dofs);
+ else
+ neighbor->get_mg_dof_indices (neighbor_dofs);
+
+ hang.distribute_local_to_global(cell_m, dofs, m);
+ // TBA Distribute neighbor dofs
+ // hang.distribute_local_to_global(cell_m, neighbor_dofs, dofs, m);
+ //
+
+ } // The neighbor is coarser
+ } // else : if face has children
+ } // skipped if no flags for face terms
+ } // faces loop
+ } // Skipped if no flags on either faces or boundary
+#else
+ //1D face terms
+#endif
+ } // cells loop
+ }
+
+ /** MG assembly routine for matrices and vectors. Hanging node
+ * constraints are taken into account. However face terms are NOT
+ * correctly handled, that is, this method is currently not
+ * capable of handling fluxex.
+
+ Different levels can be assembled. If the level is -1, then the
+ active level is used.
+ */
+ template <int dim, typename MATRIX, typename VEC>
+ void assemble(MGDoFHandler<dim> &mg_dof_handler, ConstraintMatrix &hang,
+ int level, MATRIX& m, VEC& rhs, LocalAssembleBase<dim, MGDoFHandler<dim> > &local,
+ const unsigned int this_mpi_process = 0)
+ {
+ // Cell iterators.
+ typename MGDoFHandler<dim>::active_cell_iterator cell, endc;
+ if(level < 0) {
+ cell = mg_dof_handler.begin_active();
+ endc = mg_dof_handler.end();
+ } else {
+ cell = mg_dof_handler.begin_active(level);
+ endc = mg_dof_handler.end(level);
+ }
+
+ unsigned int dofs_per_cell =
+ mg_dof_handler.get_fe().dofs_per_cell;
+
+ // dofs and neighbor dofs.
+ vector<unsigned int> dofs(dofs_per_cell);
+ vector<unsigned int> neighbor_dofs(dofs_per_cell);
+
+ FullMatrix<double> cell_m(dofs_per_cell, dofs_per_cell);
+ FullMatrix<double> cell_neighbor_m(dofs_per_cell, dofs_per_cell);
+ Vector<double> local_rhs(dofs_per_cell);
+
+ for (; cell!=endc; ++cell)
+ if( cell->subdomain_id() == this_mpi_process )
+ {
+ // Reset local matrices.
+ cell_m = 0;
+ cell_neighbor_m = 0;
+
+ // Update indices
+ if(level<0)
+ cell->get_dof_indices(dofs);
+ else
+ cell->get_mg_dof_indices (dofs);
+
+ // Assemble cell term if needed
+ if( local.flags & assemble_cell )
+ local.assemble_cell_term(cell, cell_m);
+ if( local.flags & assemble_rhs_cell)
+ local.assemble_rhs_term(cell, local_rhs);
+
+ hang.distribute_local_to_global(cell_m, dofs, m);
+ hang.distribute_local_to_global(local_rhs, dofs, rhs);
+
+#if (DEAL_II_DIM > 1)
+ // Assemble boundary or face terms if needed
+ if( local.flags & (assemble_boundary|assemble_face) ) {
+
+ for (unsigned int face_no=0;
+ face_no<GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
+ {
+ typename MGDoFHandler<dim>::face_iterator face=cell->face(face_no);
+
+ if (face->at_boundary())
+ {
+ if(local.flags & assemble_boundary) {
+ cell_m = 0;
+ local.assemble_boundary_term(cell, face_no, cell_m);
+ hang.distribute_local_to_global(cell_m, dofs, m);
+ }
+ if(local.flags & assemble_rhs_boundary) {
+ local_rhs =0;
+ local.assemble_rhs_boundary_term(cell, face_no, local_rhs);
+ hang.distribute_local_to_global(local_rhs, dofs, rhs);
+ }
+ } // End of Boundary terms
+ else if(local.flags & assemble_face) {
+ typename MGDoFHandler<dim>::cell_iterator neighbor=
+ cell->neighbor(face_no);;
+
+ if (face->has_children()) // The neighbor is more refined
+ {
+ const unsigned int n_face_no=
+ cell->neighbor_of_neighbor(face_no);
+
+
+ for (unsigned int subface_no=0;
+ subface_no<face->n_children(); ++subface_no)
+ {
+ typename MGDoFHandler<dim>::active_cell_iterator
+ neighbor_child
+ = cell->neighbor_child_on_subface (face_no, subface_no);
+
+ Assert (neighbor_child->face(n_face_no) == face->child(subface_no),
+ ExcInternalError());
+ Assert (!neighbor_child->has_children(), ExcInternalError());
+
+ cell_m = 0;
+ cell_neighbor_m = 0;
+ local.assemble_face_term(cell, face_no, subface_no,
+ neighbor_child, n_face_no,
+ cell_m, cell_neighbor_m);
+ if(level<0)
+ neighbor_child->get_dof_indices (neighbor_dofs);
+ else
+ neighbor_child->get_mg_dof_indices (neighbor_dofs);
+
+ hang.distribute_local_to_global(cell_m, dofs, m);
+ // TBA Distribute neighbors!!!
+ // hang.distribute_local_to_global(cell_m, neighbor_dofs, dofs, m);
+ //
+
+
+ }
+ } // face has children
+ else
+ {
+ if (neighbor->level() == cell->level())
+ {
+ const unsigned int n_face_no=cell->neighbor_of_neighbor(face_no);
+
+ cell_neighbor_m = 0;
+ cell_m = 0;
+ local.assemble_face_term(cell, face_no,
+ neighbor, n_face_no,
+ cell_m,
+ cell_neighbor_m);
+
+ if(level<0)
+ neighbor->get_dof_indices (neighbor_dofs);
+ else
+ neighbor->get_mg_dof_indices (neighbor_dofs);
+
+ hang.distribute_local_to_global(cell_m, dofs, m);
+ // TBA: distribute neighbor dofs....
+ // hang.distribute_local_to_global(cell_m, neighbor_dofs, dofs, m);
+ //
+
+ } // Same level on the two neighbors
+ else
+ {
+ Assert(neighbor->level() < cell->level(), ExcInternalError());
+
+ const std::pair<unsigned int, unsigned int> faceno_subfaceno=
+ cell->neighbor_of_coarser_neighbor(face_no);
+ const unsigned int neighbor_face_no=faceno_subfaceno.first,
+ neighbor_subface_no=faceno_subfaceno.second;
+
+ Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
+ neighbor_subface_no)
+ == cell,
+ ExcInternalError());
+
+ cell_neighbor_m = 0;
+ cell_m = 0;
+ local.assemble_face_term(cell, face_no,
+ neighbor, neighbor_face_no,
+ neighbor_subface_no,
+ cell_m, cell_neighbor_m);
+
+ if(level<0)
+ neighbor->get_dof_indices (neighbor_dofs);
+ else
+ neighbor->get_mg_dof_indices (neighbor_dofs);
+
+ hang.distribute_local_to_global(cell_m, dofs, m);
+ // TBA Distribute neighbor dofs
+ // hang.distribute_local_to_global(cell_m, neighbor_dofs, dofs, m);
+ //
+
+ } // The neighbor is coarser
+ } // else : if face has children
+ } // skipped if no flags for face terms
+ } // faces loop
+
+ } // Skipped if no flags on either faces or boundary
+
+#else
+ //1D face terms
+#endif
+ } // cells loop
+ }
+
+
+ /** MG assembly routine for rhs vectors. Hanging node
+ * constraints are taken into account. However face terms are NOT
+ * correctly handled, that is, this method is currently not
+ * capable of handling fluxex.
+
+ Different levels can be assembled. If the level is -1, then the
+ active level is used.
+ */
+ template <int dim, typename VECTOR>
+ void assemble_rhs(MGDoFHandler<dim> &mg_dof_handler, ConstraintMatrix &hang,
+ int level, VECTOR& rhs, LocalAssembleBase<dim, MGDoFHandler<dim> > &local,
+ const unsigned int this_mpi_process = 0)
+ {
+ // Cell iterators.
+ typename MGDoFHandler<dim>::active_cell_iterator cell, endc;
+ if(level < 0) {
+ cell = mg_dof_handler.begin_active();
+ endc = mg_dof_handler.end();
+ } else {
+ cell = mg_dof_handler.begin_active(level);
+ endc = mg_dof_handler.end(level);
+ }
+
+ unsigned int dofs_per_cell =
+ mg_dof_handler.get_fe().dofs_per_cell;
+
+ // dofs and neighbor dofs.
+ vector<unsigned int> dofs(dofs_per_cell);
+
+ Vector<double> cell_rhs(dofs_per_cell);
+
+ for (; cell!=endc; ++cell)
+ if( cell->subdomain_id() == this_mpi_process )
+ {
+ // Update indices
+ if(level<0)
+ cell->get_dof_indices(dofs);
+ else
+ cell->get_mg_dof_indices (dofs);
+
+ // Assemble cell term if needed
+ if( local.flags & assemble_rhs_cell ) {
+ // Reset the cell term
+ cell_rhs = 0;
+ local.assemble_rhs_term(cell, cell_rhs);
+ // Distribute the local rhs to the global matrix.
+ hang.distribute_local_to_global(cell_rhs, dofs, rhs);
+ }
+
+ // Assemble boundary or face terms if needed
+ if( local.flags & (assemble_rhs_boundary/* |assemble_rhs_face */) )
+ for (unsigned int face_no=0;
+ face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ {
+ typename MGDoFHandler<dim>::face_iterator face=cell->face(face_no);
+ if (face->at_boundary())
+ {
+ cell_rhs = 0;
+ local.assemble_rhs_boundary_term(cell, face_no, cell_rhs);
+ // Distribute the local rhs to the global matrix.
+ hang.distribute_local_to_global(cell_rhs, dofs, rhs);
+ }
+ } // face loop
+ } // Cell loop
+ }
+
+ /** DG-Like assembly routine for matrices and vectors. Hanging
+ * node constraints are not taken into account, and one should do
+ * this outside. However face terms are correctly handled, and
+ * this method is capable of assemblying fully DG objects. */
+ template <int dim, typename DH, typename MATRIX, typename VEC>
+ void assemble(DH &dof_handler, MATRIX& system_matrix, VEC& system_rhs,
+ LocalAssembleBase<dim, DH> &local,
+ const unsigned int this_mpi_process = 0)
+ {
+ const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+ std::vector<unsigned int> dofs (dofs_per_cell);
+ std::vector<unsigned int> dofs_neighbor (dofs_per_cell);
+
+ // Now we create the cell matrices and vectors. Here we need two
+ // cell matrices, both for face terms that include test functions
+ // <code>vi</code> (internal shape functions, i.e. shape functions
+ // of the current cell).
+ //
+ // i stands for internal, e stands for external
+ FullMatrix<double> ui_vi_matrix (dofs_per_cell, dofs_per_cell);
+ FullMatrix<double> ui_ve_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> local_rhs(dofs_per_cell);
+
+ // Furthermore we need some cell iterators.
+ typename DH::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ // Now we start the loop over all
+ // active cells.
+ for (;cell!=endc; ++cell)
+ if( cell->subdomain_id() == this_mpi_process )
+ {
+ // On each cell we need to reset the <code>ui_vi_matrix</code>
+ // and <code>cell_vector</code> to zero, before assembling the
+ // cell terms.
+ ui_vi_matrix = 0;
+ local_rhs = 0;
+
+ // and call the function that assembles the cell terms. The
+ // first argument is the <code>FEValues</code> that was
+ // previously reinit'ed on the current cell.
+ if( local.flags & assemble_cell )
+ local.assemble_cell_term(cell, ui_vi_matrix);
+
+ if( local.flags & assemble_rhs_cell )
+ local.assemble_rhs_term(cell, local_rhs);
+
+ // As in previous examples the vector `dofs' includes the
+ // dof_indices.
+ cell->get_dof_indices (dofs);
+
+ // This is the start of the nested loop over all faces.
+ if( local.flags & (assemble_boundary|assemble_face) ) {
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ {
+ // First we set the face iterator
+ typename DH::face_iterator face=cell->face(face_no);
+
+ // and clear the <code>ui_ve_matrix</code> on each face.
+ ui_ve_matrix = 0;
+
+ // Now we distinguish the four different cases in the
+ // ordering mentioned above. We start with faces belonging
+ // to the boundary of the domain.
+ if (face->at_boundary())
+ {
+ // and assemble the corresponding face terms.
+ if(local.flags & assemble_boundary)
+ local.assemble_boundary_term(cell, face_no,
+ ui_vi_matrix);
+ if(local.flags & assemble_rhs_boundary)
+ local.assemble_rhs_boundary_term(cell, face_no,
+ local_rhs);
+ }
+ else if(local.flags & assemble_face)
+ {
+ // Now we are not on the boundary of the domain,
+ // therefore there must exist a neighboring cell.
+ typename DH::cell_iterator neighbor=
+ cell->neighbor(face_no);
+
+ // We proceed with the second and most complicated case:
+ // the neighboring cell is more refined than the current
+ // cell. As in deal.II neighboring cells are restricted
+ // to have a level difference of not more than one, the
+ // neighboring cell is known to be at most ONCE more
+ // refined than the current cell. Furthermore also the
+ // face is more refined, i.e. it has children. Here we
+ // note that the following part of code will not work
+ // for <code>dim==1</code>.
+ if (face->has_children())
+ {
+ // First we store which number the current cell has
+ // in the list of neighbors of the neighboring
+ // cell. Hence,
+ // neighbor-@>neighbor(neighbor_face_no) equals the
+ // current cell <code>cell</code>.
+ const unsigned int neighbor_face_no=
+ cell->neighbor_of_neighbor(face_no);
+
+
+ // We loop over subfaces
+ for (unsigned int subface_no=0;
+ subface_no<face->n_children(); ++subface_no)
+ {
+ // and set the cell iterator
+ // <code>neighbor_child</code> to the cell
+ // placed `behind' the current subface.
+ typename DH::active_cell_iterator
+ neighbor_child
+ = cell->neighbor_child_on_subface (face_no, subface_no);
+
+ // As these are quite complicated indirections
+ // which one does not usually get right at first
+ // attempt we check for the internal
+ // consistency.
+ Assert (neighbor_child->face(neighbor_face_no) == face->child(subface_no),
+ ExcInternalError());
+ Assert (!neighbor_child->has_children(), ExcInternalError());
+
+ // We need to reset the
+ // <code>ui_ve_matrix</code> on each subface
+ // because on each subface the <code>un</code>
+ // belong to different neighboring cells.
+ ui_ve_matrix = 0;
+
+ // As already mentioned above for the current
+ // case (case 2) we employ the
+ // <code>FESubfaceValues</code> of the current
+ // cell (here reinited for the current cell,
+ // face and subface) and we employ the
+ // FEFaceValues of the neighboring child cell.
+ local.assemble_face_term(cell, face_no, subface_no,
+ neighbor_child, neighbor_face_no,
+ ui_vi_matrix,
+ ui_ve_matrix);
+
+ // Then we get the dof indices of the
+ // neighbor_child cell
+ neighbor_child->get_dof_indices (dofs_neighbor);
+
+ // and distribute <code>ui_ve_matrix</code> to
+ // the system_matrix
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ system_matrix.add(dofs_neighbor[i], dofs[k],
+ ui_ve_matrix(i,k));
+ }
+ // End of <code>if (face-@>has_children())</code>
+ }
+ else
+ {
+ // We proceed with case 3, i.e. neighboring cell is
+ // of the same refinement level as the current cell.
+ if (neighbor->level() == cell->level())
+ {
+ // Like before we store which number the current
+ // cell has in the list of neighbors of the
+ // neighboring cell.
+ const unsigned int neighbor_face_no=cell->neighbor_of_neighbor(face_no);
+
+ local.assemble_face_term(cell, face_no,
+ neighbor, neighbor_face_no,
+ ui_vi_matrix,
+ ui_ve_matrix);
+ // End of <code>if (neighbor-@>level() ==
+ // cell-@>level())</code>
+ }
+ else
+ {
+ // Finally we consider case 4. When the
+ // neighboring cell is not finer and not of the
+ // same refinement level as the current cell it
+ // must be coarser.
+ Assert(neighbor->level() < cell->level(), ExcInternalError());
+
+ // Find out the how many'th face_no and
+ // subface_no the current face is w.r.t. the
+ // neighboring cell.
+ const std::pair<unsigned int, unsigned int> faceno_subfaceno=
+ cell->neighbor_of_coarser_neighbor(face_no);
+ const unsigned int neighbor_face_no=faceno_subfaceno.first,
+ neighbor_subface_no=faceno_subfaceno.second;
+
+ Assert (neighbor->neighbor_child_on_subface (neighbor_face_no,
+ neighbor_subface_no)
+ == cell,
+ ExcInternalError());
+
+ local.assemble_face_term(cell, face_no,
+ neighbor,
+ neighbor_face_no,
+ neighbor_subface_no,
+ ui_vi_matrix,
+ ui_ve_matrix);
+ }
+
+ // Now we get the dof indices of the
+ // <code>neighbor</code> cell,
+ neighbor->get_dof_indices (dofs_neighbor);
+
+ // and distribute the
+ // <code>ui_ve_matrix</code>.
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ system_matrix.add(dofs_neighbor[i], dofs[k],
+ ui_ve_matrix(i,k));
+ }
+ // End of <code>face not at boundary</code>:
+ }
+ // End of loop over all faces:
+ }
+ }
+
+ // Finally we distribute the
+ // <code>ui_vi_matrix</code>
+ for (unsigned int i=0; i<dofs_per_cell; ++i) {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i,j));
+ system_rhs(dofs[i]) += local_rhs(i);
+ }
+ }
+ }
+
+ /** Assembly routine for rhs. Hanging node constraints are not
+ * taken into account, and one should do this outside. */
+ template <int dim, typename DH, typename VEC>
+ void assemble_rhs(DH &dof_handler, VEC& system_rhs,
+ LocalAssembleBase<dim, DH> &local,
+ const unsigned int this_mpi_process = 0)
+ {
+ const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+ std::vector<unsigned int> dofs (dofs_per_cell);
+
+ Vector<double> local_rhs(dofs_per_cell);
+
+ // Furthermore we need some cell iterators.
+ typename DH::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ // Now we start the loop over all
+ // active cells.
+ for (;cell!=endc; ++cell)
+ if( cell->subdomain_id() == this_mpi_process )
+ {
+ // On each cell we need to reset the <code>ui_vi_matrix</code>
+ // and <code>cell_vector</code> to zero, before assembling the
+ // cell terms.
+ local_rhs = 0;
+
+ // and call the function that assembles the cell terms. The
+ // first argument is the <code>FEValues</code> that was
+ // previously reinit'ed on the current cell.
+ if( local.flags & assemble_rhs_cell )
+ local.assemble_rhs_term(cell, local_rhs);
+
+ // As in previous examples the vector `dofs' includes the
+ // dof_indices.
+ cell->get_dof_indices (dofs);
+
+ // This is the start of the nested loop over all faces.
+ if( local.flags & assemble_rhs_boundary ) {
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ {
+ // First we set the face iterator
+ typename DH::face_iterator face=cell->face(face_no);
+
+ // Now we distinguish the four different cases in the
+ // ordering mentioned above. We start with faces belonging
+ // to the boundary of the domain.
+ if (face->at_boundary())
+ {
+ if(local.flags)
+ local.assemble_rhs_boundary_term(cell, face_no,
+ local_rhs);
+ }
+ }
+ }
+
+ // Finally we distribute the
+ // <code>rhs</code>
+ for (unsigned int i=0; i<dofs_per_cell; ++i) {
+ system_rhs(dofs[i]) += local_rhs(i);
+ }
+ }
+ }
+}
+
+#endif
--- /dev/null
+#ifndef LH_OUTPUT_PROCESSOR_H
+#define LH_OUTPUT_PROCESSOR_H
+
+#include <fstream>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_constraints.h>
+#include <lac/block_vector.h>
+
+#include <grid/grid_refinement.h>
+#include <grid/tria.h>
+#include <grid/filtered_iterator.h>
+
+// #include <numerics/error_estimator.h>
+#include <base/function.h>
+#include <numerics/solution_transfer.h>
+#include <numerics/data_out.h>
+
+#include <base/convergence_table.h>
+#include <base/logstream.h>
+#include <base/config.h>
+
+#include <base/parameter_handler.h>
+
+#include "vector_space.h"
+#include <map>
+
+template<int dim, typename DH=DoFHandler<dim> >
+class FilteredDataOut : public DataOut<dim, DH>
+{
+public:
+ FilteredDataOut (const unsigned int subdomain_id)
+ :
+ subdomain_id (subdomain_id)
+ {}
+
+ virtual typename DH::cell_iterator
+ first_cell ()
+ {
+ typename DH::active_cell_iterator
+ cell = this->dofs->begin_active();
+ while ((cell != this->dofs->end()) &&
+ (cell->subdomain_id() != subdomain_id))
+ ++cell;
+
+ return cell;
+ }
+
+ virtual typename DH::cell_iterator
+ next_cell (const typename DH::cell_iterator &old_cell)
+ {
+ if (old_cell != this->dofs->end())
+ {
+ const IteratorFilters::SubdomainEqualTo
+ predicate(subdomain_id);
+
+ return
+ ++(FilteredIterator
+ <typename DH::active_cell_iterator>
+ (predicate,old_cell));
+ }
+ else
+ return old_cell;
+ }
+
+private:
+ const unsigned int subdomain_id;
+};
+
+
+
+template <int dim, typename VECTOR=BlockVector<double> >
+class OutputProcessor : public Subscriptor
+{
+public:
+ /** The constructor takes the mpi initialization stuff. */
+ OutputProcessor (const unsigned int n_mpi_processes = 1,
+ const unsigned int this_mpi_process = 0);
+
+ /** Initialize the given values for the paramter file. */
+ static void declare_parameters(ParameterHandler &prm);
+
+ /** Parse the given parameter handler. */
+ void parse_parameters(ParameterHandler &prm);
+
+ /** Calculate the error of the numeric solution in variuous norms. Store
+ the result in the given table. */
+ void error_from_exact(const VectorSpace<dim> & vspace,
+ const VECTOR &solution,
+ const Function<dim> &exact,
+ unsigned int table_no = 0,
+ double dt=0.);
+
+ /** Difference between two solutions in two different vector spaces. */
+ void difference(const VectorSpace<dim> &, const VECTOR &,
+ const VectorSpace<dim> &, const VECTOR &,
+ unsigned int table_no = 0, double dt=0.);
+
+ /** Difference between two solutions in the same vector space. */
+ void difference(const VectorSpace<dim> &, const VECTOR &,
+ const VECTOR &, unsigned int table_no = 0, double dt=0.);
+
+ /** Output solution. Using some standard format. This is now deprecated. */
+ void output_solution(const VectorSpace<dim> &, const VECTOR &, const std::string &);
+
+ /** Prepare to output data on the given file. This will initialize
+ the data_out object and a file with a filename that is the
+ combination of the @p filename, eventually a
+ processor number and the output suffix. */
+ void prepare_data_output(const DoFHandler<dim> &dh, const std::string &filename);
+
+ /** Add the given vector to the output file. Prior to calling this
+ method, you have to call the prepare_data_output method. The
+ string can be a comma separated list of components, or a single
+ description. In this latter case, a progressive number per
+ component is added in the end. */
+ void add_data_vector(const VECTOR &data_vector, const std::string &desc);
+
+ /** Actually write the file. Once the data_out has been prepared,
+ vectors have been added, the data can be written to a file. This
+ is done in this class. At the end of this function call,
+ data_out and output_file are in a pristine situation, and the
+ process can be started again.*/
+ void write_data_and_clear();
+
+ /** Dump vector. Just dump the vector to a file. This is useful to
+ read back in the future and restart from where we left off.*/
+ void dump_vector(const VECTOR &, const std::string &);
+
+ /** By default output first table. */
+ void output_table(const unsigned int table_no=0);
+
+private:
+ /** Error results.*/
+ std::vector<ConvergenceTable> tables;
+
+ /** Headers for tables and output. Contains the name of the solution
+ components. */
+ std::vector<std::string> headers;
+
+ /** Headers for latex tables. Contains the name of the solution
+ components. */
+ std::vector<std::string> latex_headers;
+
+ /** Captions for latex. */
+ std::vector<std::string> latex_captions;
+
+ /** Names of the tables. */
+ std::vector<std::string> names;
+
+ /** The parameters have been read. */
+ bool initialized;
+
+ /** Write the solution. */
+ bool write_solution;
+
+ /** Solution format. */
+ std::string solution_format;
+
+ /** Output the partitioning of the domain. */
+ bool output_partitioning;
+
+ /** Number of MPI processes. */
+ const unsigned int n_mpi_processes;
+
+ /** Id of MPI process. */
+ const unsigned int this_mpi_process;
+
+ /** Output file. */
+ std::ofstream output_file;
+
+ /** Outputs only the data that refers to this process. */
+ FilteredDataOut<dim, DoFHandler<dim> > data_out;
+};
+
+#endif
--- /dev/null
+#include "../include/output_processor.h"
+
+#include <base/logstream.h>
+#include <base/quadrature_lib.h>
+#include <base/utilities.h>
+
+#include <grid/grid_tools.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <fe/mapping_q.h>
+#include <fe/fe.h>
+
+template <int dim, typename VECTOR>
+OutputProcessor<dim,VECTOR>::OutputProcessor (const unsigned int n_mpi,
+ const unsigned int this_mpi) :
+ n_mpi_processes(n_mpi),
+ this_mpi_process(this_mpi),
+ data_out(this_mpi)
+{
+ initialized = false;
+}
+
+template <int dim, typename VECTOR>
+void OutputProcessor<dim,VECTOR>::declare_parameters (ParameterHandler &prm)
+{
+ prm.declare_entry ("Solution names", "u", Patterns::Anything(),
+ "This is used to output the solution to a file.");
+
+ prm.declare_entry ("Output partitioning", "false", Patterns::Bool());
+
+ prm.enter_subsection("Solution Out Format");
+ DataOut<dim>::declare_parameters(prm);
+ prm.set("Output format", "vtk");
+ prm.leave_subsection();
+}
+template <int dim, typename VECTOR>
+void OutputProcessor<dim,VECTOR>::parse_parameters (ParameterHandler &prm)
+{
+ output_partitioning = prm.get_bool ("Output partitioning");
+
+ headers = Utilities::split_string_list(prm.get ("Solution names"));
+
+ prm.enter_subsection("Solution Out Format");
+ data_out.parse_parameters(prm);
+ prm.leave_subsection();
+
+ initialized = true;
+}
+
+
+template <int dim, typename VECTOR>
+void OutputProcessor<dim,VECTOR>::prepare_data_output(const DoFHandler<dim> &dh,
+ const std::string &filename) {
+ AssertThrow(initialized, ExcNotInitialized());
+ if(data_out.default_suffix() != "") {
+
+ // If the output is needed and we have many processes, just output
+ // the one we need *in intermediate format*.
+ std::string fname = filename;
+ if(n_mpi_processes > 1) {
+ fname += ("_" + Utilities::int_to_string(this_mpi_process, 2) +
+ data_out.default_suffix(DataOutBase::deal_II_intermediate)) ;
+ } else {
+ fname += data_out.default_suffix();
+ }
+
+ deallog << "Will write on file: " << fname.c_str() << std::endl;
+ output_file.open(fname.c_str());
+ AssertThrow(output_file, ExcIO());
+ data_out.attach_dof_handler (dh);
+
+ if(n_mpi_processes > 1) {
+ // Output the partitioning
+ if(output_partitioning) {
+ std::vector<unsigned int> partition_int (dh.get_tria().n_active_cells());
+ GridTools::get_subdomain_association (dh.get_tria(), partition_int);
+ Vector<double> partitioning(partition_int.begin(),
+ partition_int.end());
+ static Vector<double> static_partitioning;
+ static_partitioning.swap(partitioning);
+ data_out.add_data_vector (static_partitioning, "partitioning");
+ }
+ }
+ }
+}
+
+
+template <int dim, typename VECTOR>
+void OutputProcessor<dim,VECTOR>::add_data_vector(const VECTOR &data_vector,
+ const std::string &desc)
+{
+ AssertThrow(initialized, ExcNotInitialized());
+ deallog.push("AddingData");
+ std::vector<std::string> dd = Utilities::split_string_list(desc);
+ if(data_out.default_suffix() != "") {
+ if (dd.size() ==1 )
+ data_out.add_data_vector (data_vector, desc);
+ else
+ data_out.add_data_vector (data_vector, dd);
+ deallog << "Added data: " << desc << std::endl;
+ }
+ deallog.pop();
+}
+
+
+template <int dim, typename VECTOR>
+void OutputProcessor<dim,VECTOR>::write_data_and_clear() {
+ AssertThrow(initialized, ExcNotInitialized());
+ AssertThrow(output_file, ExcIO());
+ deallog.push("WritingData");
+ if(data_out.default_suffix() != "") {
+ data_out.build_patches();
+
+ if(n_mpi_processes > 1) {
+ data_out.write_deal_II_intermediate(output_file);
+ } else {
+ data_out.write(output_file);
+ }
+ deallog << "Wrote output file." << std::endl;
+ data_out.clear();
+ output_file.close();
+ deallog << "Reset output." << std::endl;
+ }
+ deallog.pop();
+}
+
+template <int dim, typename VECTOR>
+void OutputProcessor<dim,VECTOR>::dump_vector (const VECTOR &,
+ const std::string & )
+{
+ Assert(false, ExcNotImplemented());
+ // Only specializations exist.
+}
--- /dev/null
+/*
+ * Immersed Boundary Problem:
+ *
+ * Header Files
+ *
+ * Author:
+ * Luca Heltai <heltai@dimat.unipv.it>
+ * =============
+ * License: GPL.
+ * =============
+ * $Id: ibm_rhs.h,v 1.15 2005/04/05 14:46:49 luca Exp $
+ *
+ */
+#ifndef PARSED_SYMMETRIC_TENSOR_FUNCTION_H
+#define PARSED_SYMMETRIC_TENSOR_FUNCTION_H
+
+#include <base/parameter_handler.h>
+#include <base/function.h>
+#include <base/function_parser.h>
+#include <base/point.h>
+#include <base/tensor.h>
+#include <base/symmetric_tensor.h>
+
+#include <lac/vector.h>
+
+using namespace dealii;
+
+template <int rank, int dim>
+class ParsedSymmetricTensorFunction
+{
+ public:
+ /** Constructor. It initializes the function parser.*/
+ ParsedSymmetricTensorFunction();
+
+ /** Declare parameters needed by this class. */
+ static void declare_parameters(ParameterHandler &prm);
+
+ /** Parse parameters needed by this class. */
+ void parse_parameters(ParameterHandler &prm);
+
+ /** Get ONE value at the given point. Evaluate the parsed function
+ at the given point and return the symmetric tensor. */
+ const SymmetricTensor<rank, dim> & operator() (const Point<dim> &p) const;
+
+ /** Get a vector of the values at a number of specified points */
+ void value_list(const std::vector< Point<dim> > &points,
+ std::vector< SymmetricTensor<rank, dim> > &values);
+
+ /** Get the tensor currently stored. */
+ inline const SymmetricTensor<rank, dim> & operator() () const {
+ return t;
+ };
+
+ /** Get the size of the equivalent triangular matrix. If rank = 2,
+ this number is equal to dim, otherwise it is equal to dim*(dim+1)/2. */
+ static unsigned int get_dim_triangular();
+
+ /** Set internal time of the function. */
+ void set_time(double);
+ private:
+ FunctionParser<dim> f;
+ mutable SymmetricTensor<rank, dim> t;
+};
+
+
+template<int dim>
+inline
+Tensor<2,dim> operator*(const SymmetricTensor<4,dim> & C,
+ const Tensor<2,dim> &gradv) {
+ Tensor<2,dim> T;
+ for(unsigned int i=0; i<dim; ++i)
+ for(unsigned int j=0; j<dim; ++j)
+ for(unsigned int k=0; k<dim; ++k)
+ for(unsigned int l=0; l<dim; ++l)
+ T[i][j] += C[i][j][k][l] * gradv[k][l];
+ return T;
+}
+
+
+template<int dim>
+inline
+double operator*(const Tensor<2,dim> &gradv,
+ const SymmetricTensor<2,dim> & C) {
+ double T = 0;
+ for(unsigned int i=0; i<dim; ++i)
+ for(unsigned int j=0; j<dim; ++j)
+ T += C[i][j] * gradv[i][j];
+ return T;
+}
+
+
+template<int dim>
+inline
+double double_contract(const SymmetricTensor<2,dim> & C,
+ const Tensor<2,dim> &gradv) {
+ double T = 0;
+ for(unsigned int i=0; i<dim; ++i)
+ for(unsigned int j=0; j<dim; ++j)
+ T += C[i][j] * gradv[i][j];
+ return T;
+}
+
+
+
+template<int dim>
+inline
+Tensor<1,dim> operator*(const SymmetricTensor<2,dim> & C,
+ const Tensor<1,dim> &gradt) {
+ Tensor<1,dim> T;
+ for(unsigned int j=0; j<dim; ++j)
+ for(unsigned int i=0; i<dim; ++i)
+ T[j] += C[j][i] * gradt[i];
+ return T;
+}
+
+
+/* template<int dim> */
+/* inline */
+/* double double_contract(const Tensor<2,dim> &gradv, */
+/* const Tensor<2,dim> & C) { */
+/* double T = 0; */
+/* for(unsigned int i=0; i<dim; ++i) */
+/* for(unsigned int j=0; j<dim; ++j) */
+/* T += C[i][j] * gradv[i][j]; */
+/* return T; */
+/* } */
+
+
+
+/* template<int dim> */
+/* inline */
+/* Tensor<2,dim> operator*(const SymmetricTensor<4,dim> & C, */
+/* const Tensor<2,dim> &gradv) { */
+/* Tensor<2,dim> T; */
+/* for(unsigned int i=0; i<dim; ++i) */
+/* for(unsigned int j=0; j<dim; ++j) */
+/* for(unsigned int k=0; k<dim; ++k) */
+/* for(unsigned int l=0; l<dim; ++l) */
+/* T[i][j] += C[i][j][k][l] * gradv[k][l]; */
+/* return T; */
+/* } */
+
+
+/* template<int dim> */
+/* inline */
+/* double operator*(const Tensor<2,dim> &gradv, */
+/* const SymmetricTensor<2,dim> & C) { */
+/* double T = 0; */
+/* for(unsigned int i=0; i<dim; ++i) */
+/* for(unsigned int j=0; j<dim; ++j) */
+/* T += C[i][j] * gradv[i][j]; */
+/* return T; */
+/* } */
+
+
+/* template<int dim> */
+/* inline */
+/* double double_contract(const SymmetricTensor<2,dim> & C, */
+/* const Tensor<2,dim> &gradv) { */
+/* double T = 0; */
+/* for(unsigned int i=0; i<dim; ++i) */
+/* for(unsigned int j=0; j<dim; ++j) */
+/* T += C[i][j] * gradv[i][j]; */
+/* return T; */
+/* } */
+
+
+
+/* template<int dim> */
+/* inline */
+/* Tensor<1,dim> operator*(const SymmetricTensor<2,dim> & C, */
+/* const Tensor<1,dim> &gradt) { */
+/* Tensor<1,dim> T; */
+/* for(unsigned int j=0; j<dim; ++j) */
+/* for(unsigned int i=0; i<dim; ++i) */
+/* T[j] += C[j][i] * gradt[i]; */
+/* return T; */
+/* } */
+
+
+template<int dim>
+inline
+double double_contract(const Tensor<2,dim> &gradv,
+ const Tensor<2,dim> & C) {
+ double T = 0;
+ for(unsigned int i=0; i<dim; ++i)
+ for(unsigned int j=0; j<dim; ++j)
+ T += C[i][j] * gradv[i][j];
+ return T;
+}
+
+
+#endif
--- /dev/null
+/*
+ * Immersed Boundary Problem:
+ *
+ * Header Files
+ *
+ * Author:
+ * Luca Heltai <heltai@dimat.unipv.it>
+ * =============
+ * License: GPL.
+ * =============
+ * $Id: ibm_rhs.cc,v 1.34 2005/04/05 14:46:49 luca Exp $
+ *
+ */
+#include "../include/parsed_symmetric_tensor_function.h"
+#include <base/utilities.h>
+
+template <int rank, int dim>
+ void ParsedSymmetricTensorFunction<rank, dim>::value_list(const std::vector< Point<dim> > &points,
+ std::vector< SymmetricTensor<rank, dim> > &values)
+{
+
+ for (unsigned int n=0; n<points.size(); ++n)
+ values[n] = this->operator()(points[n]);
+
+}
+
+template <int rank, int dim>
+ unsigned int ParsedSymmetricTensorFunction<rank, dim>::get_dim_triangular()
+{
+ if (rank == 2)
+ return dim;
+ else if (rank == 4)
+ return (dim*(dim+1)/2);
+ else {
+ AssertThrow(false, ExcInternalError());
+ return 0;
+ }
+}
+
+
+template <int rank, int dim>
+void ParsedSymmetricTensorFunction<rank, dim>::set_time (double t) {
+ f.set_time(t);
+}
+
+
+template <int rank, int dim>
+ParsedSymmetricTensorFunction<rank, dim>::ParsedSymmetricTensorFunction () :
+ // Number of vector functions.
+ f(get_dim_triangular()*(get_dim_triangular()+1)/2)
+{
+ Assert(rank == 2 || rank == 4,
+ ExcMessage("Rank has to be even, either 2 or 4."));
+}
+
+template <int rank, int dim>
+void ParsedSymmetricTensorFunction<rank, dim>::declare_parameters(ParameterHandler &prm)
+{
+ std::string vnames;
+ switch (dim) {
+ case 1:
+ vnames = "x,t";
+ break;
+ case 2:
+ vnames = "x,y,t";
+ break;
+ case 3:
+ vnames = "x,y,z,t";
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ break;
+ }
+ prm.declare_entry("Variable names", vnames, Patterns::Anything(),
+ "The name of the variables as they will be used in the function, separated by ','.");
+ // The expressions of the function
+ std::vector<std::string> expr(get_dim_triangular());
+ for(unsigned int i=0; i<get_dim_triangular(); ++i) {
+ for(unsigned int j=i; j<get_dim_triangular(); ++j)
+ expr[i] += (i==j) ? "1" : "; 0";
+ char tmp[100];
+ sprintf(tmp, "Row %d", i+1);
+ // Now we have the function expressions
+ prm.declare_entry(tmp, expr[i], Patterns::Anything(), (i == 0) ?
+ "Separate different components expressions by ';' "
+ "as ',' is used internally by the function parser."
+ : "");
+ }
+ prm.declare_entry("Function constants", "", Patterns::Anything(),
+ "Any constant used inside the functions which is not a variable name.");
+}
+
+template <int rank, int dim>
+ void ParsedSymmetricTensorFunction<rank, dim>::parse_parameters(ParameterHandler &prm)
+{
+ std::string vnames = prm.get("Variable names");
+ std::string expr;
+ for(unsigned int i=0; i < get_dim_triangular(); ++i) {
+ char tmp[100];
+ sprintf(tmp, "Row %d", i+1);
+ expr += prm.get(tmp);
+ if(i+1 < get_dim_triangular())
+ expr += "; ";
+ }
+ std::string constants_list = prm.get("Function constants");
+
+ std::vector<std::string> const_list =
+ Utilities::split_string_list(constants_list, ',');
+ std::map<std::string, double> constants;
+ for(unsigned int i = 0; i < const_list.size(); ++i) {
+ std::vector<std::string> this_c =
+ Utilities::split_string_list(const_list[i], '=');
+ AssertThrow(this_c.size() == 2, ExcMessage("Invalid format"));
+ double tmp;
+ AssertThrow( sscanf(this_c[1].c_str(), "%lf", &tmp), ExcMessage("Double number?"));
+ constants[this_c[0]] = tmp;
+ }
+
+ constants["pi"] = M_PI;
+ constants["Pi"] = M_PI;
+
+ unsigned int nn = (Utilities::split_string_list(vnames)).size();
+ switch (nn) {
+ case dim:
+ // Time independent function
+ f.initialize(vnames, expr, constants);
+ break;
+ case dim+1:
+ // Time dependent function
+ f.initialize(vnames, expr, constants, true);
+ break;
+ default:
+ AssertThrow(false, ExcMessage("Not the correct size. Check your code."));
+ }
+}
--- /dev/null
+#ifndef UTILITIES_HLT
+#define UTILITIES_HLT
+
+#include <base/utilities.h>
+#include <grid/grid_tools.h>
+#include <lac/vector_memory.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <lac/block_vector.h>
+#include <lac/block_sparse_matrix.h>
+
+using namespace dealii;
+
+template <class Matrix>
+class InverseMatrix : public Subscriptor
+{
+ public:
+ InverseMatrix (const Matrix &m);
+
+ void vmult (Vector<double> &dst,
+ const Vector<double> &src) const;
+
+ void Tvmult (Vector<double> &dst,
+ const Vector<double> &src) const;
+
+ private:
+ const SmartPointer<const Matrix> matrix;
+
+ mutable GrowingVectorMemory<> vector_memory;
+};
+
+
+class BBt : public Subscriptor
+{
+ public:
+ void reinit (const BlockSparseMatrix<double> &A);
+
+
+ void vmult (Vector<double> &dst,
+ const Vector<double> &src) const;
+
+ void Tvmult (Vector<double> &dst,
+ const Vector<double> &src) const;
+
+ private:
+ SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+
+ mutable BlockVector<double> u, v;
+
+ unsigned int dim;
+};
+
+template <typename TYPE>
+void smart_delete (SmartPointer<TYPE> &sp) {
+ if(sp) {
+ TYPE * p = sp;
+ sp = 0;
+ delete p;
+ }
+}
+
+#endif
--- /dev/null
+#ifndef VECTOR_SPACE_H
+#define VECTOR_SPACE_H
+
+#include <base/parameter_handler.h>
+#include <fe/fe_tools.h>
+#include <grid/persistent_tria.h>
+#include <multigrid/mg_dof_handler.h>
+#include <base/smartpointer.h>
+#include <fe/mapping.h>
+#include <dofs/dof_constraints.h>
+#include <numerics/solution_transfer.h>
+
+using namespace dealii;
+
+/**
+ VectorSpace object. Anything related to the definition of the
+ Hilbert space that makes up the pde problem lives in this class.
+*/
+template <int dim>
+class VectorSpace : public Subscriptor
+{
+ public:
+
+ /** Empty constructor. */
+ VectorSpace ();
+ /** Full constructor. It initializes this object by reading a
+ parameter file and setting pointers to the given
+ triangulation. */
+ VectorSpace (ParameterHandler &prm,
+ Triangulation<dim> &dd,
+ const unsigned int n_mpi_processes=1,
+ const unsigned int this_mpi_process=0);
+
+ /** Destructor. Clears the created pointers. */
+ ~VectorSpace ();
+
+ /** Reinit. The state of this object after calling reinit, is
+ the same as after the full constructor is used.*/
+ void reinit(ParameterHandler &prm,
+ Triangulation<dim> &dd,
+ const unsigned int n_mpi_processes=1,
+ const unsigned int this_mpi_process=0,
+ const std::string space_name="Vector Space Parameters");
+
+ /** Reset the internal status of the current triangulation, leaving
+ untouched all the rest. The status of this object is as if the
+ current triangulation and dof handler had just been
+ constructed. Nothing is done to the saved triangulation, nor to
+ the saved dofhandler. */
+ void reinit();
+
+ /** Parse the parameter file. */
+ void parse_parameters(ParameterHandler &prm, const std::string space_name);
+
+ /** Generate entries in the given parameter file. */
+ static void declare_parameters(ParameterHandler &prm, const std::string space_name="Vector Space Parameters");
+
+ /** Initialize the mesh. */
+ void initialize_mesh();
+
+ /** Refine the mesh globally. */
+ void refine_globally();
+
+ /** Flag the mesh according to the given error estimator and to the
+ refinement strategy defined in the paratmeters. */
+ void flag_for_refinement(Vector<float> &error);
+
+ /** Save current vector space. */
+ void save_step();
+
+ /** Save current vector space, and reset triangulation to be pristine. */
+ void save_step_and_reinit();
+
+ /** Transfer solution from old to new grid. */
+ template <typename VECTOR>
+ void transfer(VECTOR &dst, const VECTOR &src);
+
+ /** Reorder the degrees of freedom according to the
+ parameters. This is done for the level specified on the
+ argument. If level is -1, then the active dofs are reordered. */
+ void reorder(int level = -1,
+ std::vector<unsigned int> target_comps=std::vector<unsigned int>());
+
+ /** Redistribute degrees of freedom. The optional arguments groups
+ together the components. It is useful for vector valued finite elements, when one wants to sort together*/
+ void redistribute_dofs(std::vector<unsigned int>
+ target_components=std::vector<unsigned int>());
+
+ /** Compute maximum and minimum diameter of the mesh cells. */
+ void measure_mesh();
+
+ /** Interpolate Boundary Conditions. Generates the boundary
+ conditions for the problem. This will be used after the assembly
+ procedures. */
+ void interpolate_dirichlet_bc(const Function<dim> & f,
+ std::map<unsigned int, double> & bvalues);
+
+ /** Return reference to finite element.*/
+ inline FiniteElement<dim,dim> & get_fe() {
+ return *fe;
+ }
+
+ /** Return reference to current dof handler.*/
+ inline MGDoFHandler<dim,dim> & get_dh() {
+ return *dh;
+ }
+
+ /** Return reference to previous dof handler.*/
+ inline MGDoFHandler<dim,dim> & get_last_dh() {
+ return *last_dh;
+ }
+
+ /** Return reference to current triangulation..*/
+ inline Triangulation<dim,dim> & get_tria() {
+ return *tria;
+ }
+
+ /** Return reference to current triangulation..*/
+ inline Triangulation<dim> & get_last_tria() {
+ return *last_tria;
+ }
+
+ /** Return reference to coarse triangulation..*/
+ inline Triangulation<dim> & get_coarse_tria() {
+ return *coarse;
+ }
+
+ /** Return reference to hanging node constraints..*/
+ inline ConstraintMatrix & get_hang() {
+ return hang;
+ }
+
+ /** Return reference to mapping.*/
+ inline Mapping<dim> & get_mapping() {
+ return *mapping;
+ }
+
+
+ /** Return constant reference to finite element.*/
+ inline const FiniteElement<dim,dim> & get_fe() const {
+ return *fe;
+ }
+
+ /** Return constant reference to current dof handler.*/
+ inline const MGDoFHandler<dim,dim> & get_dh() const {
+ return *dh;
+ }
+
+ /** Return constant reference to previous dof handler.*/
+ inline const MGDoFHandler<dim,dim> & get_last_dh() const {
+ return *last_dh;
+ }
+
+ /** Return reference to current triangulation..*/
+ inline const Triangulation<dim,dim> & get_tria() const {
+ return *tria;
+ }
+
+ /** Return reference to coarse triangulation..*/
+ inline const Triangulation<dim,dim> & get_coarse_tria() const {
+ return *coarse;
+ }
+
+ /** Return reference to current triangulation..*/
+ inline const Triangulation<dim,dim> & get_last_tria() const {
+ return *last_tria;
+ }
+
+ /** Return reference to hanging node constraints..*/
+ inline const ConstraintMatrix & get_hang() const {
+ return hang;
+ }
+
+ /** Return reference to mapping.*/
+ inline const Mapping<dim,dim> & get_mapping() const {
+ return *mapping;
+ }
+
+ /** Number of dofs. */
+ inline unsigned int n_dofs() {
+ return dh->n_dofs();
+ }
+
+ /** Number of dofs per process. */
+ inline const std::vector<unsigned int> & n_dofs_pp() const {
+ return local_dofs_per_process;
+ }
+
+ /** Number of dofs per block. */
+ inline const std::vector<unsigned int> & n_dofs_pb() const {
+ return dofs_per_block;
+ }
+
+ /** Number of dofs in this process. */
+ inline unsigned int n_ldofs() const {
+ return n_local_dofs;
+ }
+
+ /** Number of dofs per block per process. */
+ inline const std::vector<std::vector<unsigned int> > & n_dofs_pb_pp() const {
+ return dofs_per_block_per_process;
+ }
+
+ /** Number of dofs per process per block. */
+ inline const std::vector<std::vector<unsigned int> > & n_dofs_pp_pb() const {
+ return dofs_per_process_per_block;
+ }
+
+ /** Number of blocks. */
+ inline unsigned int n_blocks() const {
+ return number_of_blocks;
+ }
+ /** Bool that checks if local refiment is enabled.*/
+ bool enable_local_refinement;
+
+ /** Initial global refinement. */
+ unsigned int global_refinement;
+
+ /** The size of this mesh. */
+ double h_max;
+
+ /** The size of this mesh. */
+ double h_min;
+
+ /** Dirichlet Boundary Indicators. */
+ std::map<char, std::vector<bool> > dirichlet_bc;
+
+ /** Neumann Boundary Indicators. */
+ std::map<char, std::vector<bool> > neumann_bc;
+
+ /** Other Boundary Indicators. */
+ std::map<char, std::vector<bool> > other_bc;
+private:
+
+ /** For internal use. Counts the dofs per block, block and what not. */
+ void count_dofs(std::vector<unsigned int> target_components);
+
+ std::string fe_name;
+ std::vector<std::string> ordering;
+ unsigned int map_type;
+
+ std::vector<std::string> d_bc;
+ std::vector<std::string> n_bc;
+ std::vector<std::string> o_bc;
+
+ /** Pointer to the finite element used. */
+ SmartPointer<FiniteElement<dim,dim> > fe;
+
+ /** Pointer to the dofhandler used */
+ SmartPointer<MGDoFHandler<dim,dim> > dh;
+
+ /** Pointer to the last dofhandler used */
+ SmartPointer<MGDoFHandler<dim,dim> > last_dh;
+
+ /** Pointer to a pristine coarse triangulation. */
+ SmartPointer<Triangulation<dim,dim> > coarse;
+
+ /** Pointer to the current triangulation used */
+ SmartPointer<Triangulation<dim,dim> > tria;
+
+ /** Pointer to the previous triangulation used. */
+ SmartPointer<Triangulation<dim,dim> > last_tria;
+
+ /** Finite Element Mapping. Various mappings are supported. If the
+ mapping parameter is 0, then cartesian mapping is used. Else Qn, with
+ n the mapping parameter. */
+ SmartPointer<Mapping<dim,dim> > mapping;
+
+ /** Constraint Matrix. */
+ ConstraintMatrix hang;
+
+ public:
+
+ /** Number of processes. */
+ unsigned int n_mpi_processes;
+
+ /** The id of this process. */
+ unsigned int this_mpi_process;
+
+ private:
+
+ /** Local size of vectors. */
+ std::vector<unsigned int> local_dofs_per_process;
+
+ /** Local size of vectors. */
+ std::vector<unsigned int> dofs_per_block;
+
+ /** Local sizes of vectors divided by blocks. */
+ std::vector<std::vector<unsigned int> > dofs_per_block_per_process;
+
+ /** Local sizes of vectors divided by blocks. */
+ std::vector<std::vector<unsigned int> > dofs_per_process_per_block;
+
+ /** Number of local dofs. */
+ unsigned int n_local_dofs;
+
+ /** Number of Blocks. */
+ unsigned int number_of_blocks;
+
+ /** Number of cells. */
+ unsigned int number_of_cells;
+
+ /** Distortion coefficient. */
+ double distortion;
+
+ /** Refinement strategy. */
+ std::string refinement_strategy;
+
+ public:
+
+ /** Bottom fraction of refinement. */
+ double bottom_fraction;
+
+ /** Top fraction of refinement. */
+ double top_fraction;
+
+ private:
+
+ /** Maximum number of allowed cells. */
+ unsigned int max_cells;
+
+ /** The wind direction, in case we order the mesh upwind. */
+ Point<dim> wind;
+};
+#endif
--- /dev/null
+#include "../include/vector_space.h"
+#include "../include/utilities.h"
+#include <base/logstream.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_in.h>
+#include <grid/grid_out.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <fe/mapping_cartesian.h>
+#include <fe/mapping_q.h>
+#include <fe/fe.h>
+#include <base/utilities.h>
+#include <dofs/dof_renumbering.h>
+#include <dofs/dof_tools.h>
+#include <numerics/solution_transfer.h>
+#include <numerics/vectors.h>
+#include <grid/grid_tools.h>
+#include <grid/grid_refinement.h>
+
+using namespace std;
+
+template <int dim>
+VectorSpace<dim>::VectorSpace() :
+ fe(0, "Vector Space FE"),
+ dh(0, "Vector Space This DH"),
+ last_dh(0, "Vector Space Last DH"),
+ coarse(0, "Vector Space Coarse Grid"),
+ tria(0, "Vector Space This Grid"),
+ last_tria(0, "Vector Space Last Grid"),
+ mapping(0, "Vetor Space Mapping")
+{}
+
+template <int dim>
+VectorSpace<dim>::VectorSpace(ParameterHandler &prm,
+ Triangulation<dim> &dd,
+ const unsigned int n_mpi,
+ const unsigned int this_mpi) :
+ fe(0, "Vector Space FE"),
+ dh(0, "Vector Space This DH"),
+ last_dh(0, "Vector Space Last DH"),
+ coarse(0, "Vector Space Coarse Grid"),
+ tria(0, "Vector Space This Grid"),
+ last_tria(0, "Vector Space Last Grid"),
+ mapping(0, "Vetor Space Mapping"),
+ n_mpi_processes(n_mpi),
+ this_mpi_process(this_mpi)
+{
+ reinit(prm, dd, n_mpi, this_mpi);
+}
+
+template <int dim>
+VectorSpace<dim>::~VectorSpace()
+{
+ // Get rid of the dof handler object.
+ smart_delete(dh);
+ smart_delete(last_dh);
+ smart_delete(fe);
+ coarse = 0;
+ smart_delete(tria);
+ smart_delete(last_tria);
+ smart_delete(mapping);
+}
+
+template <int dim>
+void VectorSpace<dim>::save_step() {
+ deallog.push("Save");
+ // Get rid of old dh and triangulation
+ last_dh->clear();
+ last_tria->clear();
+
+ deallog << "Saving triangulation and dofs..." << std::endl;
+ // Copy old triangulation (which is the one just used)
+ last_tria->copy_triangulation(*tria);
+
+ // swap back and forth, so that we have the same dofs...
+ dh.swap(last_dh);
+ redistribute_dofs();
+ dh.swap(last_dh);
+ deallog.pop();
+}
+
+
+template <int dim>
+void VectorSpace<dim>::save_step_and_reinit() {
+ deallog.push("Save");
+ deallog << "Swapping tria and dh with old_tria and old_dh." << std::endl;
+ // Swap dofhandler.
+ last_dh.swap(dh);
+ // Swap triangulation.
+ last_tria.swap(tria);
+ // Now call reinit
+ reinit();
+ deallog.pop();
+}
+
+
+template <int dim>
+void VectorSpace<dim>::refine_globally()
+{
+ deallog.push("GlobalRefinement");
+ tria->refine_global (1);
+ redistribute_dofs();
+ deallog.pop();
+}
+
+template <int dim>
+void VectorSpace<dim>::flag_for_refinement(Vector<float> &error)
+{
+ deallog.push("Refining");
+ deallog << "Strategy: " << refinement_strategy << std::endl;
+ if(refinement_strategy == "global") {
+ tria->set_all_refine_flags();
+ } else if(refinement_strategy == "fixed_number") {
+ if(max_cells > 0)
+ GridRefinement::refine_and_coarsen_fixed_number(*tria, error,
+ top_fraction, bottom_fraction, max_cells);
+ else
+ GridRefinement::refine_and_coarsen_fixed_number(*tria, error,
+ top_fraction, bottom_fraction);
+ } else if(refinement_strategy == "fixed_fraction") {
+ if(max_cells > 0)
+ GridRefinement::refine_and_coarsen_fixed_fraction(*tria, error,
+ top_fraction, bottom_fraction, max_cells);
+ else
+ GridRefinement::refine_and_coarsen_fixed_fraction(*tria, error,
+ top_fraction, bottom_fraction);
+ } else if(refinement_strategy == "optimize") {
+ GridRefinement::refine_and_coarsen_optimize(*tria, error);
+ } else {
+ Assert(false, ExcInternalError());
+ }
+ deallog.pop();
+}
+
+
+
+template <int dim>
+void VectorSpace<dim>::reorder(int level,
+ std::vector<unsigned int> target_comps) {
+
+ if(level<0) level = dh->get_tria().n_levels()-1;
+
+ // Renumber by subdomain association only if needed
+ if(n_mpi_processes > 1) {
+ deallog << "Renumbering subdomain wise." << std::endl;
+ DoFRenumbering::subdomain_wise(static_cast<DoFHandler<dim>&>(*dh));
+ }
+
+ for(unsigned int i=0; i<ordering.size(); ++i) {
+ if(ordering[i] == "cuth") {
+ deallog << "Renumbering with Cuthill McKee algorithm." << std::endl;
+ DoFRenumbering::Cuthill_McKee(static_cast<DoFHandler<dim>&>(*dh));
+ } else if(ordering[i] == "comp") {
+ deallog << "Renumbering component wise." << std::endl;
+ DoFRenumbering::component_wise(static_cast<DoFHandler<dim>&>(*dh), target_comps);
+ } else if(ordering[i] == "upwind") {
+ DoFRenumbering::downstream(static_cast<DoFHandler<dim>&>(*dh), wind);
+ } else if(ordering[i] == "none") {}
+ else {
+ AssertThrow(false, ExcMessage("This Reordering not implemented"));
+ }
+ }
+
+}
+
+template <int dim>
+void VectorSpace<dim>::measure_mesh() {
+ typename Triangulation<dim>::active_cell_iterator cell, endc;
+ endc = tria->end();
+ h_max = 0;
+ h_min = 1000;
+ for(cell = tria->begin_active(); cell != endc; ++cell) {
+ h_max = std::max(h_max, cell->diameter());
+ h_min = std::min(h_min, cell->diameter());
+ }
+ deallog << "Max diameter of a cell: " << h_max << std::endl
+ << "Min diameter of a cell: " << h_min << std::endl;
+}
+
+template <int dim>
+void VectorSpace<dim>::reinit(ParameterHandler &prm,
+ Triangulation<dim> &dd,
+ const unsigned int n_mpi,
+ const unsigned int this_mpi,
+ const std::string space_name)
+{
+ deallog.push("VECTORSPACE");
+ n_mpi_processes = n_mpi;
+ this_mpi_process = this_mpi;
+
+ smart_delete(dh);
+ smart_delete(last_dh);
+ smart_delete(fe);
+ smart_delete(mapping);
+ smart_delete(coarse);
+ smart_delete(last_tria);
+
+ // Parse Parameters
+ parse_parameters(prm, space_name);
+
+ // generate Finite element
+ FiniteElement<dim> * fe_p=0;
+
+ try {
+ fe_p = FETools::get_fe_from_name<dim>(fe_name);
+ } catch (...) {
+
+ // Did not recognize the finite element to use. Try your own
+ // one...
+ // fe_p = new FESystem<dim>(FE_BGH<dim>(2),1, FE_DGP<dim>(1), 1);
+ throw;
+ }
+ fe = fe_p;
+ deallog << "Finite Element Space: " << fe->get_name() << endl;
+
+ if(map_type == 0) {
+ mapping= new MappingCartesian<dim>();
+ } else {
+ mapping = new MappingQ<dim>(map_type);
+ }
+
+ // Store a pointer to the triangulation
+ coarse = ⅆ
+
+ // Now generate an empty triangulation.
+ last_tria = new Triangulation<dim>();
+
+ // And a copy of the pristine version of the triangulation.
+ tria = new Triangulation<dim>();
+ tria->copy_triangulation(*coarse);
+
+ // And we initialize the current triangulation.
+ initialize_mesh();
+
+ // Generate new DoFHandlers
+ dh = new MGDoFHandler<dim>(*tria);
+ last_dh = new MGDoFHandler<dim>(*last_tria);
+
+ // Now generates the boundary maps
+ std::vector<std::vector<std::string> > bcs;
+ bcs.push_back(d_bc);
+ bcs.push_back(n_bc);
+ bcs.push_back(o_bc);
+
+ std::vector<std::map<char, std::vector<bool> > > bc_maps(3);
+
+ for(unsigned int bcnumber=0; bcnumber<3; ++bcnumber) {
+ // Get some aliases for the bcs
+ std::vector<std::string> & bc = bcs[bcnumber];
+ std::map<char, std::vector<bool> > & bc_map = bc_maps[bcnumber];
+
+ if(bc.size())
+ for(unsigned int i=0; i < bc.size(); ++i) {
+ std::vector<std::string> id_and_comps =
+ Utilities::split_string_list(bc[i], ':');
+ AssertThrow(id_and_comps.size() == 2,
+ ExcMessage("Wrong Format for boundary indicator map."));
+ std::vector<int> ids =
+ Utilities::string_to_int(Utilities::split_string_list(id_and_comps[0]));
+ std::vector<int> comps =
+ Utilities::string_to_int(Utilities::split_string_list(id_and_comps[1]));
+
+ unsigned int n_c = fe->n_components();
+ std::vector<bool> filter(n_c, false);
+ // Now check that the components make sense
+ for(unsigned int i=0; i<comps.size(); ++i) {
+ AssertThrow((unsigned int) comps[i] < n_c,
+ ExcIndexRange(comps[i], 0, n_c));
+ filter[ comps[i] ] = true;
+ }
+ // And now save these components. Merge them if the map is
+ // already stored.
+ for(unsigned int i = 0; i<ids.size(); ++i) {
+ // For each of the id, save the map just generated
+ if(bc_map.find(ids[i]) != bc_map.end()) {
+ // For each of the id, add the components that appear in
+ // the map just generated.
+ Assert(bc_map[ ids[i] ].size() == (unsigned int) n_c,
+ ExcDimensionMismatch(bc_map[ ids[i] ].size(), n_c));
+ for(unsigned int j=0; j<n_c; ++j) {
+ bc_map[ ids[i] ][j] = filter[j];
+ }
+ } else {
+ bc_map[ids[i]] = filter;
+ }
+ }
+ }
+ }
+ dirichlet_bc = bc_maps[0];
+ neumann_bc = bc_maps[1];
+ other_bc = bc_maps[2];
+
+ deallog.pop();
+}
+
+
+
+template <int dim>
+void VectorSpace<dim>::reinit() {
+ // In this case, we only get rid of the current dofhandler
+ // and triangulation, assuming they have already been saved.
+ dh->clear();
+ tria->clear();
+ // And we restore the initial copy of the triangulation.
+ tria->copy_triangulation(*coarse);
+}
+
+template <int dim>
+void VectorSpace<dim>::initialize_mesh()
+{
+ deallog.push("InitializeMesh");
+ tria->refine_global (global_refinement);
+ deallog << "Active Cells: "
+ << tria->n_active_cells()
+ << endl;
+ deallog.pop();
+}
+
+template <int dim>
+void VectorSpace<dim>::redistribute_dofs(std::vector<unsigned int> target_components ) {
+ // Repartition the triangulation
+ if(n_mpi_processes > 1)
+ GridTools::partition_triangulation (n_mpi_processes, *tria);
+
+ // Distort the mesh if necessary
+ if(distortion > 0.) tria->distort_random(distortion);
+
+ // Measure it
+ measure_mesh();
+
+ // And Distribute degrees of freedom
+ dh->distribute_dofs(*fe);
+
+ // Now count what you just distributed
+ count_dofs(target_components);
+
+ reorder(-1, target_components);
+ hang.clear();
+ DoFTools::make_hanging_node_constraints (*dh, hang);
+ hang.close();
+ deallog << "Number of constrained DOFS: " << hang.n_constraints() << std::endl;
+}
+
+template <int dim>
+void VectorSpace<dim>::declare_parameters(ParameterHandler &prm, const std::string space_name)
+{
+
+ prm.enter_subsection(space_name);
+
+ prm.declare_entry ("Finite element space", "FE_Q(1)",
+ Patterns::Anything(),
+ "The finite element space to use. For vector "
+ "finite elements use the notation "
+ "FESystem[FE_Q(2)^2-FE_DGP(1)] (e.g. Navier-Stokes). ");
+
+ prm.declare_entry ("Mapping degree", "1", Patterns::Integer(),
+ "Degree of the mapping. If 0 is used, then a Cartesian mapping is assumed.");
+ prm.declare_entry ("Dof ordering", "cuth, comp", Patterns::Anything(),
+ "Ordering of the degrees of freedom: none, comp, cuth, upwind.");
+ prm.declare_entry ("Wind direction", ".01, .01, 1", Patterns::Anything(),
+ "Direction of the wind for upwind ordering of the mesh. ");
+
+ prm.declare_entry ("Dirichlet boundary map", "1:0", Patterns::Anything(),
+ "Boundary indicator, followed by semicolomn and a list"
+ " of components to which this boundary conditions apply. "
+ "More boundary indicators can be separated by semicolumn. "
+ "1:0,1,4 ; 2,4:0,2");
+ prm.declare_entry ("Neumann boundary map", "2:0", Patterns::Anything(),
+ "Boundary indicators, followed by semicolomn and a list of "
+ "components to which this boundary conditions apply. More "
+ "boundary indicators can be separated by semicolumn. "
+ "1:0,1,4 ; 2,4:0,2");
+ prm.declare_entry ("Other boundary map", "3:0", Patterns::Anything(),
+ "Boundary indicator, followed by semicolomn and a list of "
+ "components to which this boundary conditions apply. More "
+ "boundary indicators can be separated by semicolumn. "
+ "1:0,1,4 ; 2,4:0,2");
+
+ prm.enter_subsection("Grid Parameters");
+
+ prm.declare_entry ("Global refinement", "4", Patterns::Integer());
+ prm.declare_entry ("Distortion coefficient", "0", Patterns::Double(),
+ "If this number is greater than zero, the mesh is distorted"
+ " upon refinement in order to disrupt its structureness.");
+
+ prm.declare_entry("Refinement strategy",
+ "fixed_number", Patterns::Selection("fixed_number|fixed_fraction|optimize|global"),
+ "fixed_number: the Top/Bottom threshold fraction of cells are flagged for "
+ "refinement/coarsening. "
+ "fixed_fraction: the cells whose error is Top/Bottom fraction of the total "
+ "are refined/coarsened. optmized: try to reach optimal error distribution, "
+ "assuming error is divided by 4 upon refining. global: refine all cells.");
+ prm.declare_entry("Bottom fraction", ".3", Patterns::Double());
+ prm.declare_entry("Top fraction", ".3", Patterns::Double());
+ prm.declare_entry("Max number of cells", "0", Patterns::Integer(),
+ "A number of zero means no limit. ");
+ prm.leave_subsection();
+
+
+ prm.leave_subsection();
+}
+
+ template <int dim>
+void VectorSpace<dim>::parse_parameters(ParameterHandler &prm, const std::string space_name)
+{
+ prm.enter_subsection(space_name);
+
+ fe_name = prm.get ("Finite element space");
+ map_type = prm.get_integer("Mapping degree");
+ std::string all_ordering = prm.get ("Dof ordering");
+
+ d_bc = Utilities::split_string_list(prm.get ("Dirichlet boundary map"), ';');
+ n_bc = Utilities::split_string_list(prm.get ("Neumann boundary map"), ';');
+ o_bc = Utilities::split_string_list(prm.get ("Other boundary map"), ';');
+
+ std::vector<std::string> wind_str =
+ Utilities::split_string_list(prm.get ("Wind direction") );
+ for(unsigned int i=0; (i<wind_str.size()) && (i<dim); ++i)
+ sscanf(wind_str[i].c_str(), "%lf", &wind[i]);
+
+ prm.enter_subsection("Grid Parameters");
+
+ global_refinement = prm.get_integer ("Global refinement");
+ distortion = prm.get_double("Distortion coefficient");
+
+ refinement_strategy = prm.get("Refinement strategy");
+ enable_local_refinement = !(refinement_strategy == "global");
+
+ bottom_fraction = prm.get_double("Bottom fraction");
+ top_fraction = prm.get_double("Top fraction");
+ max_cells = prm.get_integer("Max number of cells");
+ prm.leave_subsection();
+
+ prm.leave_subsection();
+
+ ordering = Utilities::split_string_list(all_ordering);
+}
+
+
+template <int dim>
+void VectorSpace<dim>::interpolate_dirichlet_bc(const Function<dim> & f,
+ std::map<unsigned int, double> & bvalues)
+{
+ deallog.push("DBC");
+ std::map<char, std::vector<bool> >::iterator
+ dmap = dirichlet_bc.begin(),
+ dmapend = dirichlet_bc.end();
+
+ unsigned int last_counted = 0;
+ for(; dmap!=dmapend; ++dmap) {
+ char id = dmap->first;
+ std::vector<bool> &filter = dmap->second;
+
+ deallog << (int) id << " :";
+ const unsigned int n_components = fe->n_components();
+ for(unsigned int i=0; i < n_components; ++i) {
+ if(filter[i]) ;//deallog << i << ", ";
+ }
+ VectorTools::interpolate_boundary_values(*mapping, (const DoFHandler<dim>&) *dh,
+ // Dirichlet boundary only...
+ id, f, bvalues, filter);
+ deallog << " #: " << bvalues.size() - last_counted << std::endl;
+ last_counted = bvalues.size();
+ }
+ deallog << "Total DBC dofs: " << bvalues.size() << std::endl;
+ deallog.pop();
+}
+
+
+
+template <int dim>
+void VectorSpace<dim>::count_dofs(std::vector<unsigned int> target_components) {
+ // Count dofs per processor
+ local_dofs_per_process.resize(n_mpi_processes);
+
+ // Output dofs per processor
+ deallog << "Number of DOFS: " << dh->n_dofs()
+ << ", (by partition:";
+ for (unsigned int p=0; p<n_mpi_processes; ++p) {
+ unsigned int nl = DoFTools::
+ count_dofs_with_subdomain_association (*dh, p);
+ local_dofs_per_process[p] = nl;
+ deallog << (p==0 ? ' ' : '+') << nl;
+ }
+ deallog << ")" << std::endl;
+
+ // Count number of blocks
+ number_of_blocks = ( target_components.size() ?
+ target_components[target_components.size()-1]+1 : 1 );
+
+ deallog << "Number of blocks: " << number_of_blocks << std::endl;
+
+ if( number_of_blocks > 1 ) {
+ // Count dofs per blocks.
+ dofs_per_block.resize( number_of_blocks );
+
+ DoFTools::count_dofs_per_component (*dh, dofs_per_block, false, target_components);
+
+ while(dofs_per_block.size() > number_of_blocks) {
+ deallog << "Old Deal.II!!!!"<< std::endl;
+ dofs_per_block.pop_back();
+ }
+
+ // Output the dofs per block
+ deallog << "Dofs by Block :";
+ for (unsigned int p=0; p<number_of_blocks; ++p) {
+ deallog << (p==0 ? ' ' : '+') << dofs_per_block[p];
+ }
+ deallog << std::endl;
+
+ // Now count dofs per block per processor
+ dofs_per_block_per_process.resize(number_of_blocks,
+ std::vector<unsigned int>(n_mpi_processes, 0));
+ dofs_per_process_per_block.resize(n_mpi_processes,
+ std::vector<unsigned int>(number_of_blocks, 0));
+
+ // We now have a vector that contains all subdomain ids for each dof
+ std::vector<unsigned int> subdomain_association(dh->n_dofs());
+ DoFTools::get_subdomain_association (*dh, subdomain_association);
+
+ // Now cycle on each block to find out who belongs here
+ for(unsigned int block = 0; block < number_of_blocks; ++block) {
+
+ // Select dofs belonging to this block
+ deallog << "Block " << block << ":";
+ std::vector<bool> mask(fe->n_components(), false);
+ for(unsigned int i=0; i<mask.size(); ++i) {
+ if(target_components[i] == block) mask[i] = true;
+ deallog << ( i==0 ? ' ' : ',' ) << mask[i];
+ }
+ deallog << ": ";
+
+ std::vector<bool> this_block_dofs(dh->n_dofs());
+ DoFTools::extract_dofs(*dh, mask, this_block_dofs);
+
+ // Reset counter
+ for(unsigned int j=0; j<n_mpi_processes; ++j) {
+ dofs_per_block_per_process[block][j] = 0;
+ dofs_per_process_per_block[j][block] = 0;
+ }
+
+ // Count dofs on this block
+ for(unsigned int i=0; i< dh->n_dofs(); ++i)
+ if(this_block_dofs[i]) {
+ ++dofs_per_block_per_process[block][subdomain_association[i]];
+ ++dofs_per_process_per_block[subdomain_association[i]][block];
+ }
+ unsigned int check_pb=0;
+ unsigned int check_pp=0;
+
+ // Output the dofs per block per processor
+ for(unsigned int j=0; j<n_mpi_processes; ++j) {
+ deallog << (j==0 ? ' ' : '+') << dofs_per_block_per_process[block][j];
+ check_pb += dofs_per_block_per_process[block][j];
+ check_pp += dofs_per_process_per_block[j][block];
+ }
+ deallog << std::endl;
+ Assert(check_pb == dofs_per_block[block],
+ ExcDimensionMismatch(check_pb, dofs_per_block[block]));
+ Assert(check_pp == dofs_per_block[block],
+ ExcDimensionMismatch(check_pp, dofs_per_block[block]));
+ }
+ } else {
+ //
+ dofs_per_block.resize(1);
+ dofs_per_block[0] = dh->n_dofs();
+ //
+ dofs_per_block_per_process.resize(1);
+ dofs_per_block_per_process[0] = local_dofs_per_process;
+ dofs_per_process_per_block.resize(n_mpi_processes);
+ for(unsigned int j=0; j<n_mpi_processes; ++j)
+ dofs_per_process_per_block[j] =
+ std::vector<unsigned int>(1, local_dofs_per_process[j]);
+
+ }
+
+ // Save the number of cells.
+ number_of_cells = tria->n_active_cells();
+ deallog << "Number of CELLS: " << number_of_cells << endl;
+
+ // Now save the number of local dofs...
+ n_local_dofs = local_dofs_per_process[this_mpi_process];
+}
--- /dev/null
+# This file was generated by the deal.II library.
+# Date = 2006/1/6
+# Time = 18: 2:34
+#
+# For a description of the UCD format see the AVS Developer's guide.
+#
+4 5 0 0 0
+1 0 0 0
+2 1 0 0
+3 0 1 0
+4 1 1 0
+1 0 quad 1 2 4 3
+2 1 line 1 2
+3 4 line 1 3
+4 2 line 2 4
+5 3 line 3 4
--- /dev/null
+#include "../include/base.templates.h"
+
+template class Base<deal_II_dimension>;
--- /dev/null
+#include "../include/camclay.templates.h"
+
+template class CamClay<deal_II_dimension>;
--- /dev/null
+#include "../include/domain.templates.h"
+
+template class Domain<deal_II_dimension>;
--- /dev/null
+#include "../include/error_handler.templates.h"
+#include <lac/block_vector.h>
+#include <lac/petsc_vector.h>
+#include <lac/vector.h>
+
+template class ErrorHandler<deal_II_dimension, Vector<double> >;
+template class ErrorHandler<deal_II_dimension, BlockVector<double> >;
+#ifdef CONTRIB_USE_PETSC
+template class ErrorHandler<deal_II_dimension, PETScWrappers::Vector >;
+#endif
+//template class ErrorHandler<deal_II_dimension, PETScWrappers::BlockVector >;
--- /dev/null
+#include "../include/linear_elastic.templates.h"
+
+
+template class LinearElastic<deal_II_dimension>;
--- /dev/null
+#include "../include/local_assemble.templates.h"
+#include <dofs/dof_handler.h>
+#include <multigrid/mg_dof_handler.h>
+
+template class LocalAssemble<deal_II_dimension, DoFHandler<deal_II_dimension> >;
+template class LocalAssemble<deal_II_dimension, MGDoFHandler<deal_II_dimension> >;
--- /dev/null
+#include "../include/local_assemble_base.templates.h"
+
+template class LocalAssembleBase<deal_II_dimension, DoFHandler<deal_II_dimension> >;
+template class LocalAssembleBase<deal_II_dimension, MGDoFHandler<deal_II_dimension> >;
+
--- /dev/null
+#include "../include/local_assemble_elastic_matrix.templates.h"
+
+
+template class LocalAssembleElasticMatrix<deal_II_dimension>;
--- /dev/null
+#include "../include/local_assemble_elastic_rhs.templates.h"
+
+
+template class LocalAssembleElasticRHS<deal_II_dimension>;
--- /dev/null
+#include "../include/local_assemble_mass.templates.h"
+#include <dofs/dof_handler.h>
+#include <multigrid/mg_dof_handler.h>
+
+template class LocalAssembleMass<deal_II_dimension, DoFHandler<deal_II_dimension> >;
+template class LocalAssembleMass<deal_II_dimension, MGDoFHandler<deal_II_dimension> >;
--- /dev/null
+#include "../include/local_assemble_plastic_project.templates.h"
+
+
+template class LocalAssemblePlasticProject<deal_II_dimension>;
--- /dev/null
+#include "../include/local_assemble_scalar_project.templates.h"
+
+
+template class LocalAssembleScalarProject<deal_II_dimension>;
--- /dev/null
+#include <iostream>
+#include "../include/base.h"
+
+int main ()
+{
+ try
+ {
+ Base<deal_II_dimension> problem;
+ problem.run ();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+
+ return 0;
+}
--- /dev/null
+#include "../include/output_processor.templates.h"
+#include <lac/block_vector.h>
+#include <lac/petsc_vector.h>
+#include <lac/vector.h>
+
+template <>
+void OutputProcessor<deal_II_dimension,Vector<double> >::dump_vector (const Vector<double> &rhs,
+ const std::string &filename )
+{
+ deallog.push("Dump");
+ deallog << "Writing: " << filename << std::endl;
+ std::ofstream out ( (filename).c_str());
+ rhs.block_write(out);
+ deallog.pop();
+}
+
+template class FilteredDataOut<deal_II_dimension, DoFHandler<deal_II_dimension> >;
+
+template class OutputProcessor<deal_II_dimension, Vector<double> >;
+template class OutputProcessor<deal_II_dimension, BlockVector<double> >;
+
+#ifdef DEAL_II_USE_PETSC
+template class OutputProcessor<deal_II_dimension, PETScWrappers::Vector >;
+#endif
+
+// template class OutputProcessor<deal_II_dimension, PETScWrappers::BlockVector >;
--- /dev/null
+/*
+ * Immersed Boundary Problem:
+ *
+ * Header Files
+ *
+ * Author:
+ * Luca Heltai <heltai@dimat.unipv.it>
+ * =============
+ * License: GPL.
+ * =============
+ * $Id: ibm_rhs.cc,v 1.34 2005/04/05 14:46:49 luca Exp $
+ *
+ */
+#include "../include/parsed_symmetric_tensor_function.templates.h"
+
+#if deal_II_dimension == 1
+template <>
+const SymmetricTensor<2, 1> &
+ ParsedSymmetricTensorFunction<2, 1>::operator()(const Point<1> &p) const
+{
+ t[0][0] = f.value(p,0);
+ return t;
+}
+
+
+template <>
+const SymmetricTensor<4, 1> &
+ ParsedSymmetricTensorFunction<4, 1>::operator()(const Point<1> &p) const
+{
+ t[0][0][0][0] = f.value(p, 0);
+ return t;
+}
+
+#endif
+
+#if deal_II_dimension == 2
+template <>
+const SymmetricTensor<2, 2> &
+ ParsedSymmetricTensorFunction<2, 2>::operator()(const Point<2> &p) const
+{
+ t[0][0] = f.value(p,0);
+ t[0][1] = f.value(p,1);
+ t[1][1] = f.value(p,2);
+ return t;
+}
+
+
+template <>
+const SymmetricTensor<4, 2> &
+ ParsedSymmetricTensorFunction<4, 2>::operator()(const Point<2> &p) const
+{
+ t[0][0][0][0] = f.value(p, 0);
+
+ t[0][0][1][1] = f.value(p, 1);
+ t[1][1][0][0] = f.value(p, 1);
+
+ t[0][0][0][1] = f.value(p, 2);
+ t[0][1][0][0] = f.value(p, 2);
+
+ t[1][1][1][1] = f.value(p, 3);
+
+ t[1][1][0][1] = f.value(p, 4);
+ t[0][1][1][1] = f.value(p, 4);
+
+ t[0][1][0][1] = f.value(p, 5);
+ t[1][0][0][1] = f.value(p, 5);
+ return t;
+}
+
+#endif
+
+#if deal_II_dimension == 3
+
+template <>
+const SymmetricTensor<2, 3> &
+ ParsedSymmetricTensorFunction<2, 3>::operator()(const Point<3> &p) const
+{
+ t[0][0] = f.value(p,0);
+ t[0][1] = f.value(p,1);
+ t[0][2] = f.value(p,2);
+
+ t[1][1] = f.value(p,3);
+ t[1][2] = f.value(p,4);
+
+ t[2][2] = f.value(p,5);
+
+ return t;
+}
+
+
+template <>
+const SymmetricTensor<4, 3> &
+ ParsedSymmetricTensorFunction<4, 3>::operator()(const Point<3> &p) const
+{
+ t[0][0][0][0] = f.value(p, 0);
+
+ t[0][0][1][1] = f.value(p, 1);
+ t[1][1][0][0] = f.value(p, 1);
+
+ t[0][0][2][2] = f.value(p, 2);
+ t[2][2][0][0] = f.value(p, 2);
+
+ t[0][0][0][1] = f.value(p, 3);
+ t[0][1][0][0] = f.value(p, 3);
+
+ t[0][0][1][2] = f.value(p, 4);
+ t[1][2][0][0] = f.value(p, 4);
+
+ t[0][0][2][0] = f.value(p, 5);
+ t[2][0][0][0] = f.value(p, 5);
+
+ t[1][1][1][1] = f.value(p, 6);
+
+ t[1][1][2][2] = f.value(p, 7);
+ t[2][2][1][1] = f.value(p, 7);
+
+ t[1][1][0][1] = f.value(p, 8);
+ t[0][1][1][1] = f.value(p, 8);
+
+ t[1][1][1][2] = f.value(p, 9);
+ t[1][2][1][1] = f.value(p, 9);
+
+ t[1][1][2][0] = f.value(p, 10);
+ t[2][0][1][1] = f.value(p, 10);
+
+ t[2][2][2][2] = f.value(p, 11);
+
+ t[2][2][0][1] = f.value(p, 12);
+ t[0][1][2][2] = f.value(p, 12);
+
+ t[2][2][1][2] = f.value(p, 13);
+ t[1][2][2][2] = f.value(p, 13);
+
+ t[2][2][0][2] = f.value(p, 14);
+ t[0][2][2][2] = f.value(p, 14);
+
+ t[0][1][0][1] = f.value(p, 15);
+
+ t[0][1][1][2] = f.value(p, 16);
+ t[1][2][0][1] = f.value(p, 16);
+
+ t[0][1][0][2] = f.value(p, 17);
+ t[0][2][0][1] = f.value(p, 17);
+
+ t[1][2][1][2] = f.value(p, 18);
+
+ t[1][2][0][2] = f.value(p, 19);
+ t[0][2][1][2] = f.value(p, 19);
+
+ t[0][2][0][2] = f.value(p, 20);
+
+ return t;
+}
+
+#endif
+
+template class ParsedSymmetricTensorFunction<2,deal_II_dimension>;
+template class ParsedSymmetricTensorFunction<4,deal_II_dimension>;
--- /dev/null
+#include "../include/vector_space.templates.h"
+
+template class VectorSpace<deal_II_dimension>;