clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2
-
-
#
#
# Usually, you will not need to change anything beyond this point.
# settings
include $D/common/Make.global_options
+CXX += -I$(VTK_DIR)/include/vtk-5.6
+LIBS += -L/usr/local/vtk/lib/vtk/ -I$(VTK_DIR)/include/vtk-5.6 \
+-lvtkCharts \
+-lvtkCommon \
+-lvtkDICOMParser \
+-lvtkFiltering \
+-lvtkGenericFiltering \
+-lvtkGeovis \
+-lvtkGraphics \
+-lvtkHybrid \
+-lvtkIO \
+-lvtkImaging \
+-lvtkInfovis \
+-lvtkNetCDF \
+-lvtkRendering \
+-lvtkViews \
+-lvtkVolumeRendering \
+-lvtkWidgets \
+-lvtkalglib \
+-lvtkexoIIc \
+-lvtkexpat \
+-lvtkfreetype \
+-lvtkftgl \
+-lvtkjpeg \
+-lvtklibxml2 \
+-lvtkmetaio \
+-lvtkpng \
+-lvtkproj4 \
+-lvtksqlite \
+-lvtksys \
+-lvtktiff \
+-lvtkverdict \
+-lvtkzlib
+
+
# Since the whole project consists of only one file, we need not
# consider difficult dependencies. We only have to declare the
# Listing of Parameters
# ---------------------
set Extend solution on the -2,2 box = true
-set External refinement = 5
+set External refinement = 7
set Number of cycles = 4
set Run 2d simulation = true
-set Run 3d simulation = true
+set Run 3d simulation = false
subsection Exact solution 2d
+++ /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$
-
-# 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 = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2
-
-# 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
-
-
-################################################################
-# This example program will only work if Trilinos is installed. If this
-# is not the case, then simply redefine the main targets to do nothing
-ifneq ($(USE_CONTRIB_TRILINOS),yes)
-default run clean:
- @echo
- @echo "==========================================================="
- @echo "= This program cannot be compiled without Trilinos. Make="
- @echo "= sure you have Trilinos installed and detected during ="
- @echo "= configuration of deal.II ="
- @echo "==========================================================="
- @echo
-
-else
-#
-################################################################
-
-
-# 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/include/deal.II/*/*.h)
-
-# Since the whole project consists of only one file, we need not
-# consider difficult dependencies. We only have to declare the
-# libraries which we want to link to the object file. deal.II has two
-# libraries: one for the debug mode version of the
-# application and one for optimized mode.
-libs.g := $(lib-deal2.g)
-libs.o := $(lib-deal2.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
-
-endif # USE_CONTRIB_TRILINOS
+++ /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 <deal.II/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 <deal.II/base/function.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/base/parsed_function.h>
-#include <deal.II/base/timer.h>
-#include <deal.II/grid/grid_tools.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/numerics/data_out.h>
-
-#include <deal.II/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.size(),
- std::vector<Tensor<1,dim> >(dim));
- std::vector< Point<dim> > points(qf.size());
-
- 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.size(); ++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.size(),
- std::vector<Tensor<1,dim> >(dim));
- std::vector< std::vector< Tensor<1,dim> > > plastic_grads(qf.size(),
- std::vector<Tensor<1,dim> >(dim));
- std::vector< Vector<double> > hard_values(qf.size(), Vector<double>(dim));
- std::vector< Point<dim> > points(qf.size());
-
- //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 <deal.II/base/logstream.h>
-#include <deal.II/base/parameter_handler.h>
-#include <deal.II/base/parsed_function.h>
-#include <deal.II/base/symmetric_tensor.h>
-#include <deal.II/base/table.h>
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/lac/sparse_direct.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/solver_control.h>
-#include <deal.II/lac/solver_gmres.h>
-#include <deal.II/fe/fe.h>
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/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.size();
-
-
- 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.size();
-
- 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 <deal.II/grid/grid_refinement.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/grid_out.h>
-#include <deal.II/base/parameter_handler.h>
-#include <deal.II/base/path_search.h>
-#include <deal.II/fe/mapping.h>
-#include <deal.II/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 <deal.II/base/logstream.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_in.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/grid/grid_tools.h>
-#include <deal.II/numerics/vectors.h>
-#include <deal.II/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 <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_constraints.h>
-#include <deal.II/lac/vector.h>
-
-#include <deal.II/grid/tria.h>
-
-// #include <deal.II/numerics/error_estimator.h>
-#include <deal.II/base/function.h>
-#include <deal.II/numerics/solution_transfer.h>
-#include <deal.II/numerics/data_out.h>
-
-#include <deal.II/base/convergence_table.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/config.h>
-
-#include <deal.II/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 <deal.II/base/logstream.h>
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/utilities.h>
-
-#include <deal.II/grid/grid_tools.h>
-
-#include <deal.II/numerics/vectors.h>
-#include <deal.II/numerics/matrices.h>
-#include <deal.II/numerics/data_out.h>
-#include <deal.II/fe/mapping_q.h>
-#include <deal.II/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 <deal.II/base/logstream.h>
-#include <deal.II/base/parameter_handler.h>
-#include <deal.II/base/parsed_function.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/sparsity_pattern.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/filtered_matrix.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/lac/sparse_ilu.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/sparse_direct.h>
-#include <deal.II/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 <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe.h>
-
-#include <fstream>
-#include <iostream>
-#include <deal.II/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 <deal.II/base/logstream.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/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 <deal.II/base/parsed_function.h>
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe.h>
-
-#include <fstream>
-#include <iostream>
-#include <deal.II/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 <deal.II/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 <deal.II/base/parsed_function.h>
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe.h>
-
-#include <fstream>
-#include <iostream>
-#include <deal.II/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 <deal.II/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 <deal.II/base/parsed_function.h>
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe.h>
-
-#include <fstream>
-#include <iostream>
-#include <deal.II/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 <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/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 <deal.II/base/quadrature_lib.h>
-#include <deal.II/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 <deal.II/base/parsed_function.h>
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe.h>
-
-#include <fstream>
-#include <iostream>
-#include <deal.II/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 <deal.II/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 <deal.II/base/parsed_function.h>
-
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe.h>
-
-#include <fstream>
-#include <iostream>
-#include <deal.II/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 <deal.II/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 <deal.II/base/logstream.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/multigrid/mg_base.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/numerics/vectors.h>
-#include <deal.II/numerics/matrices.h>
-#include <deal.II/fe/fe.h>
-#include <deal.II/base/quadrature_lib.h>
-
-#include <deal.II/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 <deal.II/dofs/dof_handler.h>
-#include <deal.II/dofs/dof_constraints.h>
-#include <deal.II/lac/block_vector.h>
-
-#include <deal.II/grid/grid_refinement.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/filtered_iterator.h>
-
-// #include <deal.II/numerics/error_estimator.h>
-#include <deal.II/base/function.h>
-#include <deal.II/numerics/solution_transfer.h>
-#include <deal.II/numerics/data_out.h>
-
-#include <deal.II/base/convergence_table.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/base/config.h>
-
-#include <deal.II/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 <deal.II/base/logstream.h>
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/utilities.h>
-
-#include <deal.II/grid/grid_tools.h>
-
-#include <deal.II/numerics/vectors.h>
-#include <deal.II/numerics/matrices.h>
-#include <deal.II/numerics/data_out.h>
-#include <deal.II/fe/mapping_q.h>
-#include <deal.II/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 <deal.II/base/parameter_handler.h>
-#include <deal.II/base/function.h>
-#include <deal.II/base/function_parser.h>
-#include <deal.II/base/point.h>
-#include <deal.II/base/tensor.h>
-#include <deal.II/base/symmetric_tensor.h>
-
-#include <deal.II/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 <deal.II/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 <deal.II/base/utilities.h>
-#include <deal.II/grid/grid_tools.h>
-#include <deal.II/lac/vector_memory.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/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 <deal.II/base/parameter_handler.h>
-#include <deal.II/fe/fe_tools.h>
-#include <deal.II/grid/persistent_tria.h>
-#include <deal.II/multigrid/mg_dof_handler.h>
-#include <deal.II/base/smartpointer.h>
-#include <deal.II/fe/mapping.h>
-#include <deal.II/dofs/dof_constraints.h>
-#include <deal.II/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 <deal.II/base/logstream.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_in.h>
-#include <deal.II/grid/grid_out.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/fe/mapping_cartesian.h>
-#include <deal.II/fe/mapping_q.h>
-#include <deal.II/fe/fe.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/dofs/dof_renumbering.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/numerics/solution_transfer.h>
-#include <deal.II/numerics/vectors.h>
-#include <deal.II/grid/grid_tools.h>
-#include <deal.II/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 <deal.II/lac/block_vector.h>
-#include <deal.II/lac/petsc_vector.h>
-#include <deal.II/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 <deal.II/dofs/dof_handler.h>
-#include <deal.II/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 <deal.II/dofs/dof_handler.h>
-#include <deal.II/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 <deal.II/lac/block_vector.h>
-#include <deal.II/lac/petsc_vector.h>
-#include <deal.II/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>;
+++ /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;
-}
// been covered in previous examples
// and will thus not be further
// commented on.
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/function.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/sparse_matrix.h>
-#include <deal.II/lac/compressed_sparsity_pattern.h>
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/precondition.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_tools.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/numerics/vectors.h>
-#include <deal.II/numerics/matrices.h>
-#include <deal.II/numerics/data_out.h>
-
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <lac/solver_cg.h>
+#include <lac/solver_gmres.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+
+#include "constrained_matrix.h"
#include <fstream>
#include <iostream>
// tri-linear elements, but we will
// now use it for bi-quadratic
// elements:
-#include <deal.II/fe/fe_q.h>
+#include <fe/fe_q.h>
// We will not read the grid from a
// file as in the previous example,
// but generate it using a function
// the solution) in each step, so we
// need the following include file
// instead of <code>grid_in.h</code>:
-#include <deal.II/grid/grid_out.h>
+#include <grid/grid_out.h>
// When using locally refined grids,
// following file contains a class
// which is used to handle these
// constraints:
-#include <deal.II/lac/constraint_matrix.h>
+#include <lac/constraint_matrix.h>
// In order to refine our grids
// locally, we need a function from
// coarsening based on the error
// indicators we have computed. This
// function is defined here:
-#include <deal.II/grid/grid_refinement.h>
+#include <grid/grid_refinement.h>
// Finally, we need a simple way to
// actually compute the refinement
// often yields quite nicely adapted
// grids for a wide class of
// problems.
-#include <deal.II/numerics/error_estimator.h>
+#include <numerics/error_estimator.h>
// Finally, this is as in previous
// programs:
// constraints due to hanging
// nodes. This is done using the
// following two function calls:
- hanging_node_constraints.condense (system_matrix);
- hanging_node_constraints.condense (system_rhs);
+ // hanging_node_constraints.condense (system_matrix);
+ // hanging_node_constraints.condense (system_rhs);
+
// Using them, degrees of freedom
// associated to hanging nodes have
// been removed from the linear
template <int dim>
void Step6<dim>::solve ()
{
- SolverControl solver_control (1000, 1e-12);
- SolverCG<> solver (solver_control);
+ SolverControl solver_control (10000, 1e-12);
+ SolverGMRES<> solver (solver_control);
- PreconditionSSOR<> preconditioner;
- preconditioner.initialize(system_matrix, 1.2);
+ PreconditionIdentity preconditioner;
+//preconditioner.initialize(system_matrix, 1.2);rGM
- solver.solve (system_matrix, solution, system_rhs,
- preconditioner);
+ double scaling = 10000;
+ typename DoFHandler<dim>::active_cell_iterator
+ cell=dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for(; cell != endc; ++cell)
+ scaling = std::min(scaling, cell->diameter()*cell->diameter());
+
+ ConstrainedOperator<Vector<double>, SparseMatrix<double> >
+ cc(system_matrix, hanging_node_constraints, scaling);
- hanging_node_constraints.distribute (solution);
+ // cc.distribute_rhs(system_rhs);
+ solution = system_rhs;
+ cc.distribute_rhs(system_rhs);
+
+ solver.solve (cc, solution, system_rhs,
+ preconditioner);
}
GridOut grid_out;
grid_out.write_eps (triangulation, output);
+
+ DataOut<dim> data_out;
+
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector (solution, "solution");
+ data_out.build_patches ();
+
+ filename = "solution-";
+ filename += ('0' + cycle);
+ filename += ".vtu";
+ output.close();
+ output.open(filename.c_str());
+ data_out.write_vtu (output);
}
data_out.add_data_vector (solution, "solution");
data_out.build_patches ();
- std::ofstream output ("final-solution.eps");
- data_out.write_eps (output);
+ std::ofstream output ("final-solution.vtk");
+ data_out.write_vtk (output);
}