From 2ac22f02a7fecc20cd4292fd518e0e2611d02a26 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Tue, 31 Aug 2010 14:24:46 +0000 Subject: [PATCH] Skeleton for quasistatic plasticity example: step-41 git-svn-id: https://svn.dealii.org/trunk@21809 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-41/CamClay.prm | 306 ++++++ deal.II/examples/step-41/Makefile | 136 +++ deal.II/examples/step-41/doc/builds-on | 1 + deal.II/examples/step-41/doc/intro.dox | 4 + deal.II/examples/step-41/doc/kind | 1 + deal.II/examples/step-41/doc/results.dox | 3 + deal.II/examples/step-41/doc/tooltip | 1 + .../examples/step-41/include/assemble_flags.h | 122 +++ deal.II/examples/step-41/include/base.h | 113 ++ .../examples/step-41/include/base.templates.h | 468 +++++++++ deal.II/examples/step-41/include/camclay.h | 172 +++ .../step-41/include/camclay.templates.h | 670 ++++++++++++ deal.II/examples/step-41/include/domain.h | 72 ++ .../step-41/include/domain.templates.h | 120 +++ .../examples/step-41/include/error_handler.h | 177 ++++ .../step-41/include/error_handler.templates.h | 396 +++++++ .../examples/step-41/include/linear_elastic.h | 77 ++ .../include/linear_elastic.templates.h | 160 +++ .../examples/step-41/include/local_assemble.h | 66 ++ .../include/local_assemble.templates.h | 67 ++ .../step-41/include/local_assemble_base.h | 102 ++ .../include/local_assemble_base.templates.h | 68 ++ .../include/local_assemble_elastic_matrix.h | 50 + .../local_assemble_elastic_matrix.templates.h | 85 ++ .../include/local_assemble_elastic_rhs.h | 65 ++ .../local_assemble_elastic_rhs.templates.h | 133 +++ .../local_assemble_hardening_project.h | 53 + .../step-41/include/local_assemble_mass.h | 45 + .../include/local_assemble_mass.templates.h | 93 ++ .../include/local_assemble_plastic_project.h | 53 + ...local_assemble_plastic_project.templates.h | 93 ++ .../include/local_assemble_scalar_project.h | 51 + .../local_assemble_scalar_project.templates.h | 84 ++ deal.II/examples/step-41/include/my_tools.h | 989 ++++++++++++++++++ .../step-41/include/output_processor.h | 176 ++++ .../include/output_processor.templates.h | 135 +++ .../parsed_symmetric_tensor_function.h | 193 ++++ ...rsed_symmetric_tensor_function.templates.h | 135 +++ deal.II/examples/step-41/include/utilities.h | 62 ++ .../examples/step-41/include/vector_space.h | 325 ++++++ .../step-41/include/vector_space.templates.h | 583 +++++++++++ deal.II/examples/step-41/mesh/square.inp | 16 + deal.II/examples/step-41/source/base.cc | 3 + deal.II/examples/step-41/source/camclay.cc | 3 + deal.II/examples/step-41/source/domain.cc | 3 + .../examples/step-41/source/error_handler.cc | 11 + .../examples/step-41/source/linear_elastic.cc | 4 + .../examples/step-41/source/local_assemble.cc | 6 + .../step-41/source/local_assemble_base.cc | 5 + .../source/local_assemble_elastic_matrix.cc | 4 + .../source/local_assemble_elastic_rhs.cc | 4 + .../step-41/source/local_assemble_mass.cc | 6 + .../source/local_assemble_plastic_project.cc | 4 + .../source/local_assemble_scalar_project.cc | 4 + deal.II/examples/step-41/source/main.cc | 37 + .../step-41/source/output_processor.cc | 26 + .../parsed_symmetric_tensor_function.cc | 158 +++ .../examples/step-41/source/vector_space.cc | 3 + 58 files changed, 7002 insertions(+) create mode 100644 deal.II/examples/step-41/CamClay.prm create mode 100644 deal.II/examples/step-41/Makefile create mode 100644 deal.II/examples/step-41/doc/builds-on create mode 100644 deal.II/examples/step-41/doc/intro.dox create mode 100644 deal.II/examples/step-41/doc/kind create mode 100644 deal.II/examples/step-41/doc/results.dox create mode 100644 deal.II/examples/step-41/doc/tooltip create mode 100644 deal.II/examples/step-41/include/assemble_flags.h create mode 100644 deal.II/examples/step-41/include/base.h create mode 100644 deal.II/examples/step-41/include/base.templates.h create mode 100644 deal.II/examples/step-41/include/camclay.h create mode 100644 deal.II/examples/step-41/include/camclay.templates.h create mode 100644 deal.II/examples/step-41/include/domain.h create mode 100644 deal.II/examples/step-41/include/domain.templates.h create mode 100644 deal.II/examples/step-41/include/error_handler.h create mode 100644 deal.II/examples/step-41/include/error_handler.templates.h create mode 100644 deal.II/examples/step-41/include/linear_elastic.h create mode 100644 deal.II/examples/step-41/include/linear_elastic.templates.h create mode 100644 deal.II/examples/step-41/include/local_assemble.h create mode 100644 deal.II/examples/step-41/include/local_assemble.templates.h create mode 100644 deal.II/examples/step-41/include/local_assemble_base.h create mode 100644 deal.II/examples/step-41/include/local_assemble_base.templates.h create mode 100644 deal.II/examples/step-41/include/local_assemble_elastic_matrix.h create mode 100644 deal.II/examples/step-41/include/local_assemble_elastic_matrix.templates.h create mode 100644 deal.II/examples/step-41/include/local_assemble_elastic_rhs.h create mode 100644 deal.II/examples/step-41/include/local_assemble_elastic_rhs.templates.h create mode 100644 deal.II/examples/step-41/include/local_assemble_hardening_project.h create mode 100644 deal.II/examples/step-41/include/local_assemble_mass.h create mode 100644 deal.II/examples/step-41/include/local_assemble_mass.templates.h create mode 100644 deal.II/examples/step-41/include/local_assemble_plastic_project.h create mode 100644 deal.II/examples/step-41/include/local_assemble_plastic_project.templates.h create mode 100644 deal.II/examples/step-41/include/local_assemble_scalar_project.h create mode 100644 deal.II/examples/step-41/include/local_assemble_scalar_project.templates.h create mode 100644 deal.II/examples/step-41/include/my_tools.h create mode 100644 deal.II/examples/step-41/include/output_processor.h create mode 100644 deal.II/examples/step-41/include/output_processor.templates.h create mode 100644 deal.II/examples/step-41/include/parsed_symmetric_tensor_function.h create mode 100644 deal.II/examples/step-41/include/parsed_symmetric_tensor_function.templates.h create mode 100644 deal.II/examples/step-41/include/utilities.h create mode 100644 deal.II/examples/step-41/include/vector_space.h create mode 100644 deal.II/examples/step-41/include/vector_space.templates.h create mode 100644 deal.II/examples/step-41/mesh/square.inp create mode 100644 deal.II/examples/step-41/source/base.cc create mode 100644 deal.II/examples/step-41/source/camclay.cc create mode 100644 deal.II/examples/step-41/source/domain.cc create mode 100644 deal.II/examples/step-41/source/error_handler.cc create mode 100644 deal.II/examples/step-41/source/linear_elastic.cc create mode 100644 deal.II/examples/step-41/source/local_assemble.cc create mode 100644 deal.II/examples/step-41/source/local_assemble_base.cc create mode 100644 deal.II/examples/step-41/source/local_assemble_elastic_matrix.cc create mode 100644 deal.II/examples/step-41/source/local_assemble_elastic_rhs.cc create mode 100644 deal.II/examples/step-41/source/local_assemble_mass.cc create mode 100644 deal.II/examples/step-41/source/local_assemble_plastic_project.cc create mode 100644 deal.II/examples/step-41/source/local_assemble_scalar_project.cc create mode 100644 deal.II/examples/step-41/source/main.cc create mode 100644 deal.II/examples/step-41/source/output_processor.cc create mode 100644 deal.II/examples/step-41/source/parsed_symmetric_tensor_function.cc create mode 100644 deal.II/examples/step-41/source/vector_space.cc diff --git a/deal.II/examples/step-41/CamClay.prm b/deal.II/examples/step-41/CamClay.prm new file mode 100644 index 0000000000..404fbd435c --- /dev/null +++ b/deal.II/examples/step-41/CamClay.prm @@ -0,0 +1,306 @@ +# 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 + + diff --git a/deal.II/examples/step-41/Makefile b/deal.II/examples/step-41/Makefile new file mode 100644 index 0000000000..7fb1643e79 --- /dev/null +++ b/deal.II/examples/step-41/Makefile @@ -0,0 +1,136 @@ +# $Id: Makefile.large 18849 2009-05-15 17:59:46Z bangerth $ + +# The large projects Makefile looks much like the one for small +# projects. Basically, only the following seven parameters need to be +# set by you: + +application-name = camclay + +# The first denotes the dimension for which the program is to be +# compiled: +deal_II_dimension = 2 + +# The second tells us the name of the executable. It is prefixed by +# `lib/' to designate its destination directory. Note that the program +# name depends on the dimension, so you can keep copies for the +# different dimensions around: +target = lib/$(application-name)-$(deal_II_dimension)d + +# The `debug-mode' variable works as in the small projects Makefile: +debug-mode = on + +# And so does the following variable. You will have to set it to +# something more reasonable, of course. +clean-up-files = + +# Finally, here is a variable which tells the `run' rule which +# parameters to pass to the executable. Usually, this will be the name +# of an input file. +run-parameters = parameter-file.prm + +# Now, this is the last variable you need to set, namely the path to +# the deal.II toplevel directory: +D = ../../ + + + +# +# +# Usually, you will not need to change anything beyond this point. +# +# +# This tells `make' where to find the global settings and rules: +include $D/common/Make.global_options + + +# First get a list of files belonging to the project. Include files +# are expected in `include/', while implementation files are expected +# in `source/'. Object files are placed into `lib/[123]d', using the +# same base name as the `.cc' file. +cc-files = $(shell echo source/*.cc) +o-files = $(cc-files:source/%.cc=lib/$(deal_II_dimension)d/%.$(OBJEXT)) +go-files = $(cc-files:source/%.cc=lib/$(deal_II_dimension)d/%.g.$(OBJEXT)) +h-files = $(wildcard include/*.h) +lib-h-files = $(shell echo $D/base/include/base/*.h \ + $D/lac/include/lac/*.h \ + $D/deal.II/include/*/*.h) + +# As before, define a list of libraries. This, of course depends on +# the dimension in which we are working: +libs.g = $(lib-deal2-$(deal_II_dimension)d.g) \ + $(lib-lac.g) \ + $(lib-base.g) +libs.o = $(lib-deal2-$(deal_II_dimension)d.o) \ + $(lib-lac.o) \ + $(lib-base.o) + + + +# Now use the information from above to define the set of libraries to +# link with and the flags to be passed to the compiler: +ifeq ($(debug-mode),on) + libraries = $(go-files) $(libs.g) + flags = $(CXXFLAGS.g) +else + libraries = $(o-files) $(libs.o) + flags = $(CXXFLAGS.o) +endif + + +# Then augment the compiler flags by a specification of the dimension +# for which the program shall be compiled: +flags += -Ddeal_II_dimension=$(deal_II_dimension) + + +# The following two rules define how to compile C++ files into object +# files: +lib/$(deal_II_dimension)d/%.g.$(OBJEXT) : + @echo =====$(application-name)=======$(deal_II_dimension)d====debug=====$(MT)== $( $@ + +include lib/Makefile.dep + diff --git a/deal.II/examples/step-41/doc/builds-on b/deal.II/examples/step-41/doc/builds-on new file mode 100644 index 0000000000..8b13789179 --- /dev/null +++ b/deal.II/examples/step-41/doc/builds-on @@ -0,0 +1 @@ + diff --git a/deal.II/examples/step-41/doc/intro.dox b/deal.II/examples/step-41/doc/intro.dox new file mode 100644 index 0000000000..446bf8b910 --- /dev/null +++ b/deal.II/examples/step-41/doc/intro.dox @@ -0,0 +1,4 @@ + +

Introduction

+ +Introduction about Cam Clay plasticity models. diff --git a/deal.II/examples/step-41/doc/kind b/deal.II/examples/step-41/doc/kind new file mode 100644 index 0000000000..56e049c91a --- /dev/null +++ b/deal.II/examples/step-41/doc/kind @@ -0,0 +1 @@ +solids diff --git a/deal.II/examples/step-41/doc/results.dox b/deal.II/examples/step-41/doc/results.dox new file mode 100644 index 0000000000..6121fb67ed --- /dev/null +++ b/deal.II/examples/step-41/doc/results.dox @@ -0,0 +1,3 @@ +

Results

+ +Astonishing results here. diff --git a/deal.II/examples/step-41/doc/tooltip b/deal.II/examples/step-41/doc/tooltip new file mode 100644 index 0000000000..03e9264e3d --- /dev/null +++ b/deal.II/examples/step-41/doc/tooltip @@ -0,0 +1 @@ +Solving Cam Clay plasticity model, getting the Nobel Prize. diff --git a/deal.II/examples/step-41/include/assemble_flags.h b/deal.II/examples/step-41/include/assemble_flags.h new file mode 100644 index 0000000000..272ddb7763 --- /dev/null +++ b/deal.II/examples/step-41/include/assemble_flags.h @@ -0,0 +1,122 @@ +//--------------------------------------------------------------------------- +// $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 +/** + * 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 + * operator | 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 ( + static_cast (f1) | + static_cast (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 + * operator & 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 ( + static_cast (f1) & + static_cast (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 diff --git a/deal.II/examples/step-41/include/base.h b/deal.II/examples/step-41/include/base.h new file mode 100644 index 0000000000..05ad294454 --- /dev/null +++ b/deal.II/examples/step-41/include/base.h @@ -0,0 +1,113 @@ +#ifndef BASE_H +#define BASE_H + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#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 +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 &my_elas, + const Vector &my_plas, + const Vector &my_other); + + void write_stress_strain (const unsigned int &cc, + const unsigned int &step, + const Vector &my_elas); + + void write_plot_values(const unsigned int &cc, + const unsigned int &step, + const Vector &my_elas, + const Vector &my_hardening, + const Vector &my_plastic, + const ParsedSymmetricTensorFunction<4,dim> &C); + + //A parameter handler object + ParameterHandler prm; + + //ErrorHandler Object + ErrorHandler > error_handler; + + //This holds the triangulation information + Domain domain; + + //A vectorspace for the thermoelasticity + VectorSpace vspace; + + //Elasticity Class - to change models, just change the class + LinearElastic elastic; + //HypoElastic elastic; + //HyperElastic elastic; + //LinearThermoElastic elastic; + + //Plasticity Class - to change models, just change the class name + CamClay plastic; + //CamClayExplicit plastic; + //CamClayInvariant plastic; + //CamClayShearBand plastic; + //J2Flow plastic; + //J2FlowThermal 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 exact_solution; + + +}; + + +#endif diff --git a/deal.II/examples/step-41/include/base.templates.h b/deal.II/examples/step-41/include/base.templates.h new file mode 100644 index 0000000000..8664990d80 --- /dev/null +++ b/deal.II/examples/step-41/include/base.templates.h @@ -0,0 +1,468 @@ +#include "../include/base.h" + +template +Base::Base () + : + exact_solution(dim) +{} + +template +Base::~Base () +{} + +template +void Base::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 +void Base::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 +void Base::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 < +void Base::elastic_predictor (double &time) +{ + deallog.push("ELASTIC"); + + elastic.reinit_step(time); + + elastic.build_rhs(vspace); + + elastic.solve(vspace, lin_red_tol); + + deallog.pop(); +} + +template +void Base::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 +void Base::write_files (const unsigned int &cc, + const unsigned int &step, + const Vector &my_elas, + const Vector &my_plas, + const Vector &my_other) +{ + DataOut data_out; + data_out.attach_dof_handler (vspace.get_dh()); + + std::vector elastic_names; + std::vector plastic_names; + std::vector 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 partition_int (vspace.get_tria().n_active_cells()); + GridTools::get_subdomain_association (vspace.get_tria(), partition_int); + const Vector 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 +void Base::write_stress_strain (const unsigned int &cc, + const unsigned int &step, + const Vector &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 qf(QTrapez<1>(),2); + FEValues fe_v(vspace.get_fe(), qf, + update_gradients | update_quadrature_points); + std::vector< std::vector< Tensor<1,dim> > > total_grads(qf.n_quadrature_points, + std::vector >(dim)); + std::vector< Point > points(qf.n_quadrature_points); + + output << "SCALARS uyy double 1" << std::endl; + output << "LOOKUP_TABLE default" << std::endl; + + typename MGDoFHandler::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 +void Base::write_plot_values(const unsigned int &cc, + const unsigned int &step, + const Vector &my_elas, + const Vector &my_hardening, + const Vector &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(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 qf(QTrapez<1>(),1); + QTrapez qf; + FEValues fe_v(vspace.get_fe(), qf, + update_values | update_gradients | update_quadrature_points); + std::vector< std::vector< Tensor<1,dim> > > elastic_grads(qf.n_quadrature_points, + std::vector >(dim)); + std::vector< std::vector< Tensor<1,dim> > > plastic_grads(qf.n_quadrature_points, + std::vector >(dim)); + std::vector< Vector > hard_values(qf.n_quadrature_points, Vector(dim)); + std::vector< Point > points(qf.n_quadrature_points); + + //std::vector< std::vector > strain(dim, std::vector(dim) ); + + //calculate p,q,deve + typename MGDoFHandler::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 +void Base::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::declare_parameters(prm, dim); + prm.leave_subsection(); + + // ============================== + + prm.read_input("CamClay.prm"); + + // ============================== + //Initialize the domain + domain.reinit(prm); + /* + domain.get_tria().clear(); + + Point lowerleft(0.0,0.0); + Point upperright(0.3,0.1); + std::vector 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(); + +} + + diff --git a/deal.II/examples/step-41/include/camclay.h b/deal.II/examples/step-41/include/camclay.h new file mode 100644 index 0000000000..e926280a71 --- /dev/null +++ b/deal.II/examples/step-41/include/camclay.h @@ -0,0 +1,172 @@ +#ifndef CAMCLAY +#define CAMCLAY + +//deal.ii packages +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +//Trilinos packages +#include + +//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 +class CamClay +{ + public: + + CamClay(); + + ~CamClay(); + + void declare_parameters(ParameterHandler &prm); + + void parse_parameters(ParameterHandler &prm); + + void reinit(VectorSpace &vspace); + + void reinit_step(); + + void project_strain(ParameterHandler &prm, + VectorSpace &vspace, + SparseMatrix &A); + + void project_hardening(VectorSpace &vspace); + + void build_matrix(VectorSpace &vspace); + + void initial_conditions(VectorSpace &vspace); + + void update_internal_variables(VectorSpace &vspace, + Vector &elastic_solution, + ParsedSymmetricTensorFunction<4,dim> &C); + + void compute_jacobian(FullMatrix &jac, + const Vector &sol, + const double p, + const double q, + const double M, + const SymmetricTensor<4,dim> &C_qp, + const SymmetricTensor<2,dim> &xi); + + Vector dF_dstress(const SymmetricTensor<2,dim> &stress, + const double k); + + //this is the plastic strain projected into the displacement + Vector solution; + + //this is the average number of iterations per cell for the current step + Vector iterations; + + //internal variable vectors + Table<2, SymmetricTensor<2,dim> > plastic_strain; + Table<2, double > hardening; + Table<2, double> iter_table; + + Vector 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 &res, + FullMatrix &jac, + const Vector &sol, + const Vector &prev_step, + const SymmetricTensor<2,dim> &total_strain); + /* + //computes the residual + void compute_residual(Vector &res, + const Vector &prev_step, + const Vector &sol, + const SymmetricTensor<2,dim> &stress); + + + + //creates the jacobian numerical + void numerical_jacobian(FullMatrix &jac, + const Vector &sol, + const Vector &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 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 yield_stress; + + //some constants + //double shear_mod; + //double bulk_mod; + double M; + + SymmetricTensor<4,dim> C_qp; + + SparseMatrix MM; + + SparsityPattern sp_MM; + + + +}; + +#endif diff --git a/deal.II/examples/step-41/include/camclay.templates.h b/deal.II/examples/step-41/include/camclay.templates.h new file mode 100644 index 0000000000..da3c49385f --- /dev/null +++ b/deal.II/examples/step-41/include/camclay.templates.h @@ -0,0 +1,670 @@ +#include "../include/camclay.h" + +template +CamClay::CamClay() + : +yield_stress(1) +{} + +template +CamClay::~CamClay() +{ + MM.clear(); +} + +template +void CamClay::declare_parameters(ParameterHandler &prm) +{ + prm.enter_subsection("Yield Stress"); + Functions::ParsedFunction::declare_parameters(prm, 1); + prm.leave_subsection(); + +} + +template +void CamClay::parse_parameters(ParameterHandler &prm) +{ + prm.enter_subsection("Yield Stress"); + yield_stress.parse_parameters(prm); + prm.leave_subsection(); +} + +template +void CamClay::reinit(VectorSpace &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 &> (vspace.get_dh()), sp_MM); + sp_MM.compress(); + + MM.reinit(sp_MM); +} + +template +void CamClay::reinit_step() +{ + + solution = 0; + iterations = 0; + sol_hard_iter = 0; + +} + +template +void CamClay::initial_conditions(VectorSpace &vspace) +{ + + QGauss qf_v((2*vspace.get_fe().degree) + 1); + + FEValues fe_v (vspace.get_fe(), qf_v, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int n_qp_v = qf_v.n_quadrature_points; + + + typename MGDoFHandler::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; qpindex(), qp) = 2e8; + + //plastic_strain(cell->index(), qp)[0][1] = 0; + + + } + } + + + +} + +template +void CamClay::update_internal_variables(VectorSpace &vspace, + Vector &elastic_solution, + ParsedSymmetricTensorFunction<4, dim> &C) +{ + + //first thing, we need to loop over the cells + //and then the quadrature points + QGauss qf_v((2*vspace.get_fe().degree) + 1); + + FEValues fe_v (vspace.get_fe(), qf_v, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int n_qp_v = qf_v.n_quadrature_points; + + std::vector > points(n_qp_v); + + SymmetricTensor<2,dim> trial_strain; + SymmetricTensor<2,dim> trial_stress; + SymmetricTensor<2,dim> total_strain; + + std::vector< std::vector > > total_grads(n_qp_v, std::vector >(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::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 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 +inline double CamClay::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 +inline SymmetricTensor<2,3> CamClay::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 +inline double CamClay::q(const SymmetricTensor<2, 3> &xi) +{ + + double qval = std::sqrt(3.0/2.0) * xi.norm(); + + return qval; + +} + +template +double CamClay::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 cell_res(size); + Vector cell_sol(size); + Vector cell_prev(size); + FullMatrix 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= 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 +void CamClay::compute_res_jac(Vector &res, + FullMatrix &jac, + const Vector &sol, + const Vector &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 fad_double; + unsigned int size = sol.size(); + + std::vector< std::vector > estrain(dim, std::vector(dim)); + + //vectors of the unknowns and independent vars + std::vector x(size); //this is ep11,ep22,ep12,k,dgamma + std::vector off_diag(((dim*dim)-dim)/2); + std::vector 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) { + off_diag[i-size] = sol(i-3); + off_diag[i-size].diff(i,(2*(size-1))-dim); + } + + } + + for(unsigned int i=0; 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 > stress(dim, std::vector(dim)); + fad_double p; + p=0; + for(unsigned int i=0; i > s(dim, std::vector(dim)); + fad_double q; + fad_double trs; + trs = 0; + + for (unsigned int i=0; i 1e-8) deallog << "Shit is Broke" < +inline std::vector CamClay::voigt2sym(const unsigned int i) +{ + + std::vector 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 +inline double CamClay::delta(const unsigned int i, + const unsigned int j) +{ + double delta_val=0; + + if(i == j) delta_val = 1; + + return delta_val; + +} + +template +inline double CamClay::delta(const unsigned int i) +{ + double delta_val=1; + + if(i >= dim) delta_val = 0; + + return delta_val; + +} + + +template +inline unsigned int CamClay::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 +double CamClay::h(const double p, + const double k) +{ + + double hval = k * ( 2*p - k ); + + return hval; + +} + +template +Vector CamClay::dF_dstress(const SymmetricTensor<2, dim> &stress, + const double k) +{ + SymmetricTensor<2,dim> tmp; + Vector tmpv(((dim*dim)+dim)/2); + + double p = (1.0/3.0) * first_invariant(stress); + + tmp = stress; + for(unsigned int i=0; i +inline double CamClay::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 +void CamClay::project_strain(ParameterHandler &prm, + VectorSpace &vspace, + SparseMatrix &A) +{ + + LocalAssemblePlasticProject 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 +void CamClay::build_matrix(VectorSpace &vspace) +{ + QGauss quad((2*vspace.get_fe().degree) + 1); + + MatrixCreator::create_mass_matrix(vspace.get_dh(), quad, MM); + +} + +template +void CamClay::project_hardening(VectorSpace &vspace) +{ + //we need a mass matrix for a scalar variable over the mesh + //so, lets build one from scratch + + LocalAssembleScalarProject 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); + + + +} + diff --git a/deal.II/examples/step-41/include/domain.h b/deal.II/examples/step-41/include/domain.h new file mode 100644 index 0000000000..c8a2064678 --- /dev/null +++ b/deal.II/examples/step-41/include/domain.h @@ -0,0 +1,72 @@ +#ifndef DOMAIN_H +#define DOMAIN_H +#include + +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + +/** + Domain object. +*/ +template +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 & 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 tria; + /** Helper class to output the grid. */ + GridOut gridout; +}; +#endif diff --git a/deal.II/examples/step-41/include/domain.templates.h b/deal.II/examples/step-41/include/domain.templates.h new file mode 100644 index 0000000000..a1663ac830 --- /dev/null +++ b/deal.II/examples/step-41/include/domain.templates.h @@ -0,0 +1,120 @@ +#include "domain.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +template +Domain::Domain() : + search_mesh("MESH", 1) +{ + initialized = false; +} + +template +Domain::Domain(ParameterHandler &prm) : + search_mesh("MESH", 1) +{ + reinit(prm); +} + +template +Domain::~Domain() +{ + // tria.clear(); +} + + +template +void Domain::reinit(ParameterHandler &prm) +{ + deallog.push("DOMAIN"); + parse_parameters(prm); + deallog << "Generating coarse triangulation." << endl; + create_mesh(); + output_mesh(); + deallog.pop(); +} + + +template +void Domain::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::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 +void Domain::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 +void Domain::create_mesh() +{ + if(read_mesh) { + GridIn 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::parse_format(input_mesh_format)); + } else { + Point corner; + for(unsigned int d=0; d(), corner, true); +// GridTools::partition_triangulation (n_mpi_processes, tria); + } + initialized = true; +} + + +template +void Domain::output_mesh(std::ostream &out) const { + Assert(initialized, ExcNotInitialized()); + gridout.write (tria, out); +} + + +template +void Domain::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(); + } +} diff --git a/deal.II/examples/step-41/include/error_handler.h b/deal.II/examples/step-41/include/error_handler.h new file mode 100644 index 0000000000..239a251980 --- /dev/null +++ b/deal.II/examples/step-41/include/error_handler.h @@ -0,0 +1,177 @@ +#ifndef LH_ERROR_HANDLER_H +#define LH_ERROR_HANDLER_H + +#include + +#include +#include +#include + +#include + +// #include +#include +#include +#include + +#include +#include +#include + +#include + +#include + +enum NormFlags { + None = 0x00, + Linfty = 0x01, + L2 = 0x02, + W1infty = 0x04, + H1 = 0x08, + AddUp = 0x10 +}; + +using namespace dealii; + +template > +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 & vspace, + const VEC &solution, + const Function &exact, + unsigned int table_no = 0, + double dt=0.); + + /** Difference between two solutions in two different vector spaces. */ + void difference(const DoFHandler &, const VEC &, + const DoFHandler &, const VEC &, + unsigned int table_no = 0, double dt=0.); + + /** Difference between two solutions in the same vector space. */ + void difference(const DoFHandler &, 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 tables; + + /** Headers for tables and output. Contains the name of the solution + components. */ + std::vector headers; + + /** Headers for latex tables. Contains the name of the solution + components. */ + std::vector latex_headers; + + /** Captions for latex. */ + std::vector latex_captions; + + /** Names of the tables. */ + std::vector names; + + /** Type of error to compute per components. */ + std::vector > 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 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 > extras; + + /** Wether or not to calculate the rates according to the given keys. */ + std::vector 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 + * operator | 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 ( + static_cast (f1) | + static_cast (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 + * operator & 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 ( + static_cast (f1) & + static_cast (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 diff --git a/deal.II/examples/step-41/include/error_handler.templates.h b/deal.II/examples/step-41/include/error_handler.templates.h new file mode 100644 index 0000000000..93816d8567 --- /dev/null +++ b/deal.II/examples/step-41/include/error_handler.templates.h @@ -0,0 +1,396 @@ +#include "../include/error_handler.h" + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template +ErrorHandler::ErrorHandler () +{ + initialized = false; +} + +template +void ErrorHandler::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 +void ErrorHandler::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 (headers.size())); + add_rates.resize(names.size()); + tables.resize(names.size()); + latex_captions.resize(names.size()); + std::map 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 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 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_types = + Utilities::split_string_list(all_comps[j]); + for(unsigned int k=0; k +void ErrorHandler::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 +void ErrorHandler::difference(const DoFHandler & 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(0, headers.size()), table_no, dt); +} + + +template +void ErrorHandler::error_from_exact(const DoFHandler & dh, + const VECTOR &solution, + const Function &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 > error( exact.n_components, std::vector(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 select_component ( component, 1. , exact.n_components); + + Vector difference_per_cell (dh.get_tria().n_active_cells()); + + QGauss 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#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 +class LinearElastic +{ + public: + + LinearElastic(); + + ~LinearElastic(); + + void declare_parameters(ParameterHandler &prm); + + void parse_parameters(ParameterHandler &prm); + + void reinit(VectorSpace &vspace); + + void build_matrix(ParameterHandler &prm, + VectorSpace &vspace); + + void reinit_step(double &time); + + void build_rhs(VectorSpace &vspace); + + void solve(VectorSpace &vspace, + double &tolerance); + + SparseMatrix A; + + ParsedSymmetricTensorFunction<4, dim> C; + + //The solution vectors + Vector sol_total; + Vector sol_increment; + + private: + + //various functions for the laplace equation + ParsedFunction dbc; + ParsedFunction nbc; + ParsedFunction bf; + + //The system matrix for linear elasticity + SparsityPattern sp_A; + + //The distributed right hand side + Vector rhs; + + + +}; + +#endif diff --git a/deal.II/examples/step-41/include/linear_elastic.templates.h b/deal.II/examples/step-41/include/linear_elastic.templates.h new file mode 100644 index 0000000000..02b51cf708 --- /dev/null +++ b/deal.II/examples/step-41/include/linear_elastic.templates.h @@ -0,0 +1,160 @@ +#include "../include/linear_elastic.h" + +template +LinearElastic::LinearElastic() + : + dbc(dim), + nbc(dim), + bf(dim) +{} + + +template +LinearElastic::~LinearElastic() +{ + A.clear(); +} + +template +void LinearElastic::declare_parameters(ParameterHandler &prm) +{ + prm.enter_subsection("Dirichlet Data"); + Functions::ParsedFunction::declare_parameters(prm, dim); + prm.leave_subsection(); + + prm.enter_subsection("f - Body Source"); + Functions::ParsedFunction::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::declare_parameters(prm, dim); + prm.leave_subsection(); + +} + +template +void LinearElastic::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 +void LinearElastic::reinit(VectorSpace &vspace) +{ + + sp_A.reinit(vspace.n_dofs(), + vspace.n_dofs(), + vspace.get_dh().max_couplings_between_dofs()); + + DoFTools::make_sparsity_pattern (static_cast &> (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 +void LinearElastic::build_matrix(ParameterHandler &prm, + VectorSpace &vspace) +{ + + LocalAssembleElasticMatrix 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 +void LinearElastic::reinit_step(double &time) +{ + + rhs = 0; + //sol_increment = 0; + bf.set_time(time); + nbc.set_time(time); + dbc.set_time(time); +} + +template +void LinearElastic::build_rhs(VectorSpace &vspace) +{ + + LocalAssembleElasticRHS 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 +void LinearElastic::solve(VectorSpace &vspace, + double &tolerance) +{ + + //get the DBC map + std::map 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 > 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 > precon; + precon.initialize(A, 0.8); + FilteredMatrix > filtered_precon(precon); + + SolverControl control (vspace.n_dofs(), + tolerance*rhs.linfty_norm(), + false, true); + + GrowingVectorMemory > mem; + SolverCG > 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; + +} + diff --git a/deal.II/examples/step-41/include/local_assemble.h b/deal.II/examples/step-41/include/local_assemble.h new file mode 100644 index 0000000000..26f5cf631d --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble.h @@ -0,0 +1,66 @@ +#ifndef LOCAL_ASSEMBLE_STANDARD +#define LOCAL_ASSEMBLE_STANDARD + +#include "local_assemble_base.h" +#include +#include +#include +#include + +#include +#include +#include + + +template > + class LocalAssemble : public LocalAssembleBase +{ +public: +#if deal_II_dimension != 1 + LocalAssemble(); + + /** Given two face finite elements, assemble the two matrices.*/ + virtual void assemble_face_terms (FEFaceValuesBase &fe_v, + FEFaceValuesBase &fe_n_v, + FullMatrix &, + FullMatrix &); + + virtual void assemble_face_term + (const typename DH::active_cell_iterator& , + const unsigned int, + const typename DH::active_cell_iterator& , + const unsigned int, + FullMatrix &, + FullMatrix &); + + 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 &, + FullMatrix &); + + 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 &, + FullMatrix &); + + /** A pointer to fe_face_values objects. */ + SmartPointer > fe_face_v; + + /** A pointer to fe_face_values objects. */ + SmartPointer > fe_face_n_v; + + /** A pointer to fe_face_values objects. */ + SmartPointer > fe_sub_face_v; +#endif + +}; + +#endif diff --git a/deal.II/examples/step-41/include/local_assemble.templates.h b/deal.II/examples/step-41/include/local_assemble.templates.h new file mode 100644 index 0000000000..1a2d4a598e --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble.templates.h @@ -0,0 +1,67 @@ +#include "local_assemble.h" +#include + +#if deal_II_dimension != 1 + +template +LocalAssemble::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 +void LocalAssemble::assemble_face_terms (FEFaceValuesBase &, + FEFaceValuesBase &, + FullMatrix &, + FullMatrix &) +{ + Assert(false, ExcPureFunctionCalled()); +} + + +template +void LocalAssemble::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 &cell_m, + FullMatrix &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 +void LocalAssemble::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 &cell_m, + FullMatrix &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 +void LocalAssemble::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 &cell_m, + FullMatrix &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 diff --git a/deal.II/examples/step-41/include/local_assemble_base.h b/deal.II/examples/step-41/include/local_assemble_base.h new file mode 100644 index 0000000000..6e42d4de30 --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble_base.h @@ -0,0 +1,102 @@ +#ifndef __deal2__local_assemble_base_h +#define __deal2__local_assemble_base_h + +#include +#include +#include +#include +#include "assemble_flags.h" + +#include +#include + +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 > +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& fe_v, + // FullMatrix &u_v_matrix) const; + virtual void assemble_cell_term + (const typename DH::active_cell_iterator&, + FullMatrix &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 &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 &cell_m, + FullMatrix &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 &cell_m, + FullMatrix &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 &cell_m, + FullMatrix &neighbor_cell_m); + + /** Assemble rhs. This will be called for each cell.*/ + virtual void assemble_rhs_term + (const typename DH::active_cell_iterator&, + Vector &); + + /** 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 &); + + AssembleFlags flags; +}; + +#endif diff --git a/deal.II/examples/step-41/include/local_assemble_base.templates.h b/deal.II/examples/step-41/include/local_assemble_base.templates.h new file mode 100644 index 0000000000..98fdc4e873 --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble_base.templates.h @@ -0,0 +1,68 @@ +#include "local_assemble_base.h" + +template +void LocalAssembleBase::assemble_cell_term (const typename DH::active_cell_iterator&, + FullMatrix &) +{ + AssertThrow(false, ExcPureFunctionCalled()); +} + +template +void LocalAssembleBase::assemble_boundary_term +(const typename DH::active_cell_iterator&,const unsigned int, + FullMatrix &) +{ + AssertThrow(false, ExcPureFunctionCalled()); +} + +template +void LocalAssembleBase::assemble_face_term +(const typename DH::active_cell_iterator&,const unsigned int, + const typename DH::active_cell_iterator&,const unsigned int, + FullMatrix &, + FullMatrix &) +{ + AssertThrow(false, ExcPureFunctionCalled()); +} + +template +void LocalAssembleBase::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 &, + FullMatrix &) +{ + AssertThrow(false, ExcPureFunctionCalled()); +} + +template +void LocalAssembleBase::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 &, + FullMatrix &) +{ + AssertThrow(false, ExcPureFunctionCalled()); +} + + +template +void LocalAssembleBase::assemble_rhs_term +(const typename DH::active_cell_iterator&, + Vector &) +{ + AssertThrow(false, ExcPureFunctionCalled()); +} + +template +void LocalAssembleBase::assemble_rhs_boundary_term +(const typename DH::active_cell_iterator&, + const unsigned, + Vector &) +{ + AssertThrow(false, ExcPureFunctionCalled()); +} + diff --git a/deal.II/examples/step-41/include/local_assemble_elastic_matrix.h b/deal.II/examples/step-41/include/local_assemble_elastic_matrix.h new file mode 100644 index 0000000000..41c096ad27 --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble_elastic_matrix.h @@ -0,0 +1,50 @@ +#ifndef LOCAL_ASSEMBLE_ELASTIC_MATRIX +#define LOCAL_ASSEMBLE_ELASTIC_MATRIX + +#include "local_assemble.h" +#include "parsed_symmetric_tensor_function.h" +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +template +class LocalAssembleElasticMatrix : public LocalAssembleBase +{ + public: + + LocalAssembleElasticMatrix(); + + ~LocalAssembleElasticMatrix(); + + virtual void assemble_cell_term + (const typename MGDoFHandler::active_cell_iterator& cell, + FullMatrix &); + + + void reinit(FiniteElement&); + + void parameters(ParameterHandler &prm); + + private: + /** A pointer to the finite element.*/ + SmartPointer > fe; + + /** A pointer to the rhs function.*/ + ParsedSymmetricTensorFunction<4, dim> C; + + /** A pointer to fe_values objects. */ + SmartPointer > fe_v; + +}; + +#endif diff --git a/deal.II/examples/step-41/include/local_assemble_elastic_matrix.templates.h b/deal.II/examples/step-41/include/local_assemble_elastic_matrix.templates.h new file mode 100644 index 0000000000..3841fe3684 --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble_elastic_matrix.templates.h @@ -0,0 +1,85 @@ +#include "../include/local_assemble_elastic_matrix.h" +#include + +template +LocalAssembleElasticMatrix::LocalAssembleElasticMatrix() : +fe(0, "Local Assemble Fe Pointer"), + fe_v(0, "Local Assemble FeValues Pointer") +{ +} + + +template +void LocalAssembleElasticMatrix::reinit (FiniteElement& myfe) +{ + if(fe_v) { + FEValues * p = fe_v; + fe_v = 0; + delete p; + } + fe = &myfe; + QGauss quadrature(2*fe->degree + 1); + QGauss face_quadrature(2*fe->degree + 1); + UpdateFlags flags (update_values | + update_gradients | + update_q_points | + update_JxW_values); + fe_v = new FEValues(*fe, quadrature, flags); + + this->flags = assemble_cell; + +} + +template +LocalAssembleElasticMatrix::~LocalAssembleElasticMatrix() +{ + if(fe_v) { + FEValues * p = fe_v; + fe_v = 0; + delete p; + } + fe = 0; +} + +template +void LocalAssembleElasticMatrix::parameters(ParameterHandler &prm) +{ + prm.enter_subsection("Elastic Moduli"); + C.parse_parameters(prm); + prm.leave_subsection(); +} + +template +void LocalAssembleElasticMatrix::assemble_cell_term +(const typename MGDoFHandler::active_cell_iterator& cell, + FullMatrix &cell_m) +{ + cell_m = 0; + Assert(fe, ExcNotInitialized()); + + fe_v->reinit(cell); + + std::vector > points(fe_v->n_quadrature_points); + points = fe_v->get_quadrature_points(); + unsigned int comp_i, comp_j; + + for(unsigned int point =0; pointn_quadrature_points; ++point) { + SymmetricTensor<4,dim> C_qp = C(points[point]); + + for (unsigned int i=0; idofs_per_cell; ++i) { + comp_i = fe->system_to_component_index(i).first; + for (unsigned int j=0; jdofs_per_cell; ++j) { + comp_j = fe->system_to_component_index(j).first; + + for(unsigned int b=0; bshape_grad(j, point)[n] * + fe_v->shape_grad(i, point)[b] * + fe_v->JxW(point); + } //n + }//b + } //j + }//i + }//qp +} diff --git a/deal.II/examples/step-41/include/local_assemble_elastic_rhs.h b/deal.II/examples/step-41/include/local_assemble_elastic_rhs.h new file mode 100644 index 0000000000..d8bbbf3d91 --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble_elastic_rhs.h @@ -0,0 +1,65 @@ +#ifndef LOCAL_ASSEMBLE_ELASTIC_RHS +#define LOCAL_ASSEMBLE_ELASTIC_RHS + +#include "local_assemble.h" +#include "parsed_symmetric_tensor_function.h" +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +template +class LocalAssembleElasticRHS : public LocalAssembleBase +{ + public: + + LocalAssembleElasticRHS(); + + ~LocalAssembleElasticRHS(); + + virtual void assemble_rhs_term + (const typename MGDoFHandler::active_cell_iterator&, + Vector &); + + /** This object will be called for each boundary face of a triangulation.*/ + virtual void assemble_rhs_boundary_term + (const typename MGDoFHandler::active_cell_iterator&, const unsigned int, + Vector &); + + + void reinit(FiniteElement&, + Function &, + Function &, + std::map > &); + + private: + /** A pointer to the finite element.*/ + SmartPointer > fe; + + /** A pointer to the rhs function.*/ + SmartPointer > rhs; + + /** A pointer to the Neumann function.*/ + SmartPointer > neumann; + + /** A map of ids and components for neumann bc. */ + std::map > neumann_map; + + /** A pointer to fe_values objects. */ + SmartPointer > fe_v; + + /** A pointer to fe_face_values objects. */ + SmartPointer > fe_face_v; + +}; + +#endif diff --git a/deal.II/examples/step-41/include/local_assemble_elastic_rhs.templates.h b/deal.II/examples/step-41/include/local_assemble_elastic_rhs.templates.h new file mode 100644 index 0000000000..8954b90454 --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble_elastic_rhs.templates.h @@ -0,0 +1,133 @@ +#include "../include/local_assemble_elastic_rhs.h" +#include + +template +LocalAssembleElasticRHS::LocalAssembleElasticRHS() : +fe(0, "Local Assemble Fe Pointer"), + fe_v(0, "Local Assemble FeValues Pointer"), + fe_face_v(0, "Local Assemble FeFaceValues Pointer") +{ +} + + +template +void LocalAssembleElasticRHS::reinit (FiniteElement &myfe, + Function &bf, + Function &nbc, + std::map > & n_map) +{ + if(fe_v) { + FEValues * p = fe_v; + fe_v = 0; + delete p; + } + if(fe_face_v) { + FEFaceValues * p = fe_face_v; + fe_face_v = 0; + delete p; + } + fe = &myfe; + QGauss quadrature(2*fe->degree + 1); + QGauss face_quadrature(2*fe->degree + 1); + UpdateFlags flags (update_values | + update_gradients | + update_q_points | + update_JxW_values); + fe_v = new FEValues(*fe, quadrature, flags); + + fe_face_v = new FEFaceValues(*fe, face_quadrature, flags); + + this->flags = + assemble_rhs_cell| + assemble_rhs_boundary; + + rhs = &bf; + neumann = &nbc; + neumann_map = n_map; +} + +template +LocalAssembleElasticRHS::~LocalAssembleElasticRHS() +{ + if(fe_v) { + FEValues * p = fe_v; + fe_v = 0; + delete p; + } + if(fe_face_v) { + FEFaceValues * p = fe_face_v; + fe_face_v = 0; + delete p; + } + fe = 0; +} + + +template +void LocalAssembleElasticRHS::assemble_rhs_term +(const typename MGDoFHandler::active_cell_iterator& cell, + Vector &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 > load_vector (n_q_points, Vector(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; idofs_per_cell; ++i) { + comp_i = fe->system_to_component_index(i).first; + for (unsigned int q_point=0; q_pointshape_value(i,q_point) * + fe_v->JxW(q_point) ); + } + } +} + +template +void LocalAssembleElasticRHS::assemble_rhs_boundary_term +(const typename MGDoFHandler::active_cell_iterator& cell, + const unsigned int face_no, + Vector &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 & 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 > neumann_vector (n_q_points, Vector(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; idofs_per_cell; ++i) { + comp_i = fe->system_to_component_index(i).first; + if(filter[comp_i]) + for (unsigned int q_point=0; q_pointshape_value(i,q_point) * + fe_face_v->JxW(q_point) ); + } + } +} diff --git a/deal.II/examples/step-41/include/local_assemble_hardening_project.h b/deal.II/examples/step-41/include/local_assemble_hardening_project.h new file mode 100644 index 0000000000..4fa1de4237 --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble_hardening_project.h @@ -0,0 +1,53 @@ +#ifndef LOCAL_ASSEMBLE_HARDENING_PROJECT +#define LOCAL_ASSEMBLE_HARDENING_PROJECT + +#include "local_assemble.h" +#include "parsed_symmetric_tensor_function.h" +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +template +class LocalAssemblePlasticProject : public LocalAssembleBase +{ + public: + + LocalAssemblePlasticProject(); + + ~LocalAssemblePlasticProject(); + + virtual void assemble_rhs_term + (const typename MGDoFHandler::active_cell_iterator&, + Vector &); + + void parameters(ParameterHandler &prm); + + + void reinit(FiniteElement&, + Table<2, SymmetricTensor<2,dim> > &); + + private: + /** A pointer to the finite element.*/ + SmartPointer > fe; + + /** A pointer to the plastic strain */ + Table > plastic_strain; + + /** A pointer to fe_values objects. */ + SmartPointer > fe_v; + + ParsedSymmetricTensorFunction<4, dim> C; + +}; + +#endif diff --git a/deal.II/examples/step-41/include/local_assemble_mass.h b/deal.II/examples/step-41/include/local_assemble_mass.h new file mode 100644 index 0000000000..501d967777 --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble_mass.h @@ -0,0 +1,45 @@ +#ifndef LOCAL_ASSEMBLE_MASS_MATRIX +#define LOCAL_ASSEMBLE_MASS_MATRIX + +#include "local_assemble.h" + +#include +#include +#include +#include +#include +#include + +template +class LocalAssembleMass : +public LocalAssembleBase +{ + public: + + LocalAssembleMass(); + + ~LocalAssembleMass(); + + virtual void assemble_cell_term + (const typename DH::active_cell_iterator& cell, + FullMatrix &); + + virtual void assemble_rhs_term + (const typename DH::active_cell_iterator&, + Vector &); + + void reinit(FiniteElement&, + Function &); + + private: + /** A pointer to the finite element.*/ + SmartPointer > fe; + + /** A pointer to the rhs function.*/ + SmartPointer > rhs; + + /** A pointer to fe_values objects. */ + SmartPointer > fe_v; +}; + +#endif diff --git a/deal.II/examples/step-41/include/local_assemble_mass.templates.h b/deal.II/examples/step-41/include/local_assemble_mass.templates.h new file mode 100644 index 0000000000..cf33a5637a --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble_mass.templates.h @@ -0,0 +1,93 @@ +#include "../include/local_assemble_mass.h" +#include "utilities.h" +#include +#include + +template +LocalAssembleMass::LocalAssembleMass() : + fe(0, "Local Assemble Fe Pointer"), + fe_v(0, "Local Assemble FeValues Pointer") +{ +} + + +template +void LocalAssembleMass::reinit (FiniteElement &myfe, + Function &f) +{ + smart_delete(fe_v); + fe = &myfe; + QGauss quadrature(2*fe->degree + 1); + UpdateFlags flags (update_values | + update_q_points | + update_JxW_values); + fe_v = new FEValues(*fe, quadrature, flags); + + this->flags = + assemble_cell| + assemble_rhs_cell; + + rhs = &f; +} + +template +LocalAssembleMass::~LocalAssembleMass() +{ + smart_delete(fe_v); + fe = 0; +} + +template +void LocalAssembleMass::assemble_cell_term +(const typename DH::active_cell_iterator& cell, + FullMatrix &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; idofs_per_cell; ++i) { + comp_i = fe->system_to_component_index(i).first; + for (unsigned int j=0; jdofs_per_cell; ++j) { + comp_j = fe->system_to_component_index(j).first; + if(comp_i == comp_j) + for(unsigned int q_point =0; q_pointn_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 +void LocalAssembleMass::assemble_rhs_term +(const typename DH::active_cell_iterator& cell, + Vector &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 > load_vector (n_q_points, Vector(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; idofs_per_cell; ++i) { + comp_i = fe->system_to_component_index(i).first; + for (unsigned int q_point=0; q_pointshape_value(i,q_point) * + fe_v->JxW(q_point) ); + } + } +} diff --git a/deal.II/examples/step-41/include/local_assemble_plastic_project.h b/deal.II/examples/step-41/include/local_assemble_plastic_project.h new file mode 100644 index 0000000000..4a8710c847 --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble_plastic_project.h @@ -0,0 +1,53 @@ +#ifndef LOCAL_ASSEMBLE_PLASTIC_PROJECT +#define LOCAL_ASSEMBLE_PLASTIC_PROJECT + +#include "local_assemble.h" +#include "parsed_symmetric_tensor_function.h" +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +template +class LocalAssemblePlasticProject : public LocalAssembleBase +{ + public: + + LocalAssemblePlasticProject(); + + ~LocalAssemblePlasticProject(); + + virtual void assemble_rhs_term + (const typename MGDoFHandler::active_cell_iterator&, + Vector &); + + void parameters(ParameterHandler &prm); + + + void reinit(FiniteElement&, + Table<2, SymmetricTensor<2,dim> > &); + + private: + /** A pointer to the finite element.*/ + SmartPointer > fe; + + /** A pointer to the plastic strain */ + Table > plastic_strain; + + /** A pointer to fe_values objects. */ + SmartPointer > fe_v; + + ParsedSymmetricTensorFunction<4, dim> C; + +}; + +#endif diff --git a/deal.II/examples/step-41/include/local_assemble_plastic_project.templates.h b/deal.II/examples/step-41/include/local_assemble_plastic_project.templates.h new file mode 100644 index 0000000000..4b74e2d923 --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble_plastic_project.templates.h @@ -0,0 +1,93 @@ +#include "../include/local_assemble_plastic_project.h" +#include + +template +LocalAssemblePlasticProject::LocalAssemblePlasticProject() : +fe(0, "Local Assemble Fe Pointer"), + fe_v(0, "Local Assemble FeValues Pointer") +{ +} + + +template +void LocalAssemblePlasticProject::reinit (FiniteElement &myfe, + Table<2, SymmetricTensor<2,dim> > &my_plastic_strain) +{ + if(fe_v) { + FEValues * p = fe_v; + fe_v = 0; + delete p; + } + fe = &myfe; + QGauss quadrature(2*fe->degree + 1); + QGauss face_quadrature(2*fe->degree + 1); + UpdateFlags flags (update_gradients | + update_q_points | + update_JxW_values); + fe_v = new FEValues(*fe, quadrature, flags); + + this->flags = assemble_rhs_cell; + + plastic_strain = my_plastic_strain; +} + +template +LocalAssemblePlasticProject::~LocalAssemblePlasticProject() +{ + if(fe_v) { + FEValues * p = fe_v; + fe_v = 0; + delete p; + } + fe = 0; +} + +template +void LocalAssemblePlasticProject::parameters(ParameterHandler &prm) +{ + prm.enter_subsection("Elastic Moduli"); + C.parse_parameters(prm); + prm.leave_subsection(); +} + + +template +void LocalAssemblePlasticProject::assemble_rhs_term +(const typename MGDoFHandler::active_cell_iterator& cell, + Vector &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 > load_vector (n_q_points, Vector(size) ); + std::vector > points = fe_v->get_quadrature_points(); + + for (unsigned int qp=0; qp C_qp = C(points[qp]); + + for (unsigned int i=0; idofs_per_cell; ++i) { + unsigned int comp_i = fe->system_to_component_index(i).first; + + for(unsigned int b=0; bindex(), qp)[m][n] * + fe_v->shape_grad(i, qp)[b] * + fe_v->JxW(qp) ); + + } //n + } //m + } //b + } //i + } //qp +} + + diff --git a/deal.II/examples/step-41/include/local_assemble_scalar_project.h b/deal.II/examples/step-41/include/local_assemble_scalar_project.h new file mode 100644 index 0000000000..baf9515923 --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble_scalar_project.h @@ -0,0 +1,51 @@ +#ifndef LOCAL_ASSEMBLE_SCALAR_PROJECT +#define LOCAL_ASSEMBLE_SCALAR_PROJECT + +#include "local_assemble.h" +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +template +class LocalAssembleScalarProject : public LocalAssembleBase +{ + public: + + LocalAssembleScalarProject(); + + ~LocalAssembleScalarProject(); + + virtual void assemble_rhs_term + (const typename MGDoFHandler::active_cell_iterator&, + Vector &); + + + void reinit(FiniteElement&, + Table<2, double > &, + Table<2, double > &); + + private: + /** A pointer to the finite element.*/ + SmartPointer > fe; + + /** A pointer to the plastic strain */ + Table hard_table; + Table iter_table; + + /** A pointer to fe_values objects. */ + SmartPointer > fe_v; + + +}; + +#endif diff --git a/deal.II/examples/step-41/include/local_assemble_scalar_project.templates.h b/deal.II/examples/step-41/include/local_assemble_scalar_project.templates.h new file mode 100644 index 0000000000..63118f3937 --- /dev/null +++ b/deal.II/examples/step-41/include/local_assemble_scalar_project.templates.h @@ -0,0 +1,84 @@ +#include "../include/local_assemble_scalar_project.h" +#include + +template +LocalAssembleScalarProject::LocalAssembleScalarProject() : +fe(0, "Local Assemble Fe Pointer"), + fe_v(0, "Local Assemble FeValues Pointer") +{ +} + + +template +void LocalAssembleScalarProject::reinit (FiniteElement &myfe, + Table<2, double> &my_hard, + Table<2, double> &my_iter) +{ + if(fe_v) { + FEValues * p = fe_v; + fe_v = 0; + delete p; + } + fe = &myfe; + QGauss quadrature(2*fe->degree + 1); + QGauss face_quadrature(2*fe->degree + 1); + UpdateFlags flags (update_values | + update_q_points | + update_JxW_values); + fe_v = new FEValues(*fe, quadrature, flags); + + this->flags = assemble_rhs_cell; + + hard_table = my_hard; + iter_table = my_iter; +} + +template +LocalAssembleScalarProject::~LocalAssembleScalarProject() +{ + if(fe_v) { + FEValues * p = fe_v; + fe_v = 0; + delete p; + } + fe = 0; +} + + +template +void LocalAssembleScalarProject::assemble_rhs_term +(const typename MGDoFHandler::active_cell_iterator& cell, + Vector &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; qpdofs_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 +} + + diff --git a/deal.II/examples/step-41/include/my_tools.h b/deal.II/examples/step-41/include/my_tools.h new file mode 100644 index 0000000000..22796325cc --- /dev/null +++ b/deal.II/examples/step-41/include/my_tools.h @@ -0,0 +1,989 @@ +//--------------------------------------------------------------------------- +// $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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#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 + void assemble(DH &dof_handler, MATRIX& system_matrix, + LocalAssembleBase &local, + const unsigned int this_mpi_process = 0) + { + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + std::vector dofs (dofs_per_cell); + std::vector 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 + // vi (internal shape functions, i.e. shape functions + // of the current cell). + // + // i stands for internal, e stands for external + FullMatrix ui_vi_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix 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 ui_vi_matrix + // and cell_vector 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 FEValues 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::faces_per_cell; ++face_no) + { + // First we set the face iterator + typename DH::face_iterator face=cell->face(face_no); + + // and clear the ui_ve_matrix 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 dim==1. + 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 cell. + const unsigned int neighbor_face_no= + cell->neighbor_of_neighbor(face_no); + + + // We loop over subfaces + for (unsigned int subface_no=0; + subface_non_children(); ++subface_no) + { + // and set the cell iterator + // neighbor_child 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 + // ui_ve_matrix on each subface + // because on each subface the un + // belong to different neighboring cells. + ui_ve_matrix = 0; + + // As already mentioned above for the current + // case (case 2) we employ the + // FESubfaceValues 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 ui_ve_matrix to + // the system_matrix + for (unsigned int i=0; iif (face-@>has_children()) + } + 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 if (neighbor-@>level() == + // cell-@>level()) + } + 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 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 + // neighbor cell, + neighbor->get_dof_indices (dofs_neighbor); + + // and distribute the + // ui_ve_matrix. + for (unsigned int i=0; iface not at boundary: + } + // End of loop over all faces: + } + } + + // Finally we distribute the + // ui_vi_matrix + for (unsigned int i=0; i + void assemble(MGDoFHandler &mg_dof_handler, ConstraintMatrix &hang, + int level, MATRIX& m, LocalAssembleBase > &local, + const unsigned int this_mpi_process = 0) + { + AssertThrow( !(local.flags & assemble_face), ExcNotImplemented()); + // Cell iterators. + typename MGDoFHandler::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 dofs(dofs_per_cell); + vector neighbor_dofs(dofs_per_cell); + + FullMatrix cell_m(dofs_per_cell, dofs_per_cell); + FullMatrix 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 1) + // Assemble boundary or face terms if needed + if( local.flags & (assemble_boundary|assemble_face) ) { + + for (unsigned int face_no=0; + face_no::faces_per_cell; + ++face_no) + { + typename MGDoFHandler::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::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_non_children(); ++subface_no) + { + typename MGDoFHandler::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 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 + void assemble(MGDoFHandler &mg_dof_handler, ConstraintMatrix &hang, + int level, MATRIX& m, VEC& rhs, LocalAssembleBase > &local, + const unsigned int this_mpi_process = 0) + { + // Cell iterators. + typename MGDoFHandler::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 dofs(dofs_per_cell); + vector neighbor_dofs(dofs_per_cell); + + FullMatrix cell_m(dofs_per_cell, dofs_per_cell); + FullMatrix cell_neighbor_m(dofs_per_cell, dofs_per_cell); + Vector 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::faces_per_cell; + ++face_no) + { + typename MGDoFHandler::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::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_non_children(); ++subface_no) + { + typename MGDoFHandler::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 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 + void assemble_rhs(MGDoFHandler &mg_dof_handler, ConstraintMatrix &hang, + int level, VECTOR& rhs, LocalAssembleBase > &local, + const unsigned int this_mpi_process = 0) + { + // Cell iterators. + typename MGDoFHandler::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 dofs(dofs_per_cell); + + Vector 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::faces_per_cell; ++face_no) + { + typename MGDoFHandler::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 + void assemble(DH &dof_handler, MATRIX& system_matrix, VEC& system_rhs, + LocalAssembleBase &local, + const unsigned int this_mpi_process = 0) + { + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + std::vector dofs (dofs_per_cell); + std::vector 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 + // vi (internal shape functions, i.e. shape functions + // of the current cell). + // + // i stands for internal, e stands for external + FullMatrix ui_vi_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix ui_ve_matrix (dofs_per_cell, dofs_per_cell); + Vector 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 ui_vi_matrix + // and cell_vector 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 FEValues 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::faces_per_cell; ++face_no) + { + // First we set the face iterator + typename DH::face_iterator face=cell->face(face_no); + + // and clear the ui_ve_matrix 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 dim==1. + 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 cell. + const unsigned int neighbor_face_no= + cell->neighbor_of_neighbor(face_no); + + + // We loop over subfaces + for (unsigned int subface_no=0; + subface_non_children(); ++subface_no) + { + // and set the cell iterator + // neighbor_child 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 + // ui_ve_matrix on each subface + // because on each subface the un + // belong to different neighboring cells. + ui_ve_matrix = 0; + + // As already mentioned above for the current + // case (case 2) we employ the + // FESubfaceValues 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 ui_ve_matrix to + // the system_matrix + for (unsigned int i=0; iif (face-@>has_children()) + } + 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 if (neighbor-@>level() == + // cell-@>level()) + } + 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 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 + // neighbor cell, + neighbor->get_dof_indices (dofs_neighbor); + + // and distribute the + // ui_ve_matrix. + for (unsigned int i=0; iface not at boundary: + } + // End of loop over all faces: + } + } + + // Finally we distribute the + // ui_vi_matrix + for (unsigned int i=0; i + void assemble_rhs(DH &dof_handler, VEC& system_rhs, + LocalAssembleBase &local, + const unsigned int this_mpi_process = 0) + { + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + std::vector dofs (dofs_per_cell); + + Vector 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 ui_vi_matrix + // and cell_vector to zero, before assembling the + // cell terms. + local_rhs = 0; + + // and call the function that assembles the cell terms. The + // first argument is the FEValues 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::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 + // rhs + for (unsigned int i=0; i + +#include +#include +#include + +#include +#include +#include + +// #include +#include +#include +#include + +#include +#include +#include + +#include + +#include "vector_space.h" +#include + +template > +class FilteredDataOut : public DataOut +{ +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 + + (predicate,old_cell)); + } + else + return old_cell; + } + +private: + const unsigned int subdomain_id; +}; + + + +template > +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 & vspace, + const VECTOR &solution, + const Function &exact, + unsigned int table_no = 0, + double dt=0.); + + /** Difference between two solutions in two different vector spaces. */ + void difference(const VectorSpace &, const VECTOR &, + const VectorSpace &, const VECTOR &, + unsigned int table_no = 0, double dt=0.); + + /** Difference between two solutions in the same vector space. */ + void difference(const VectorSpace &, 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 &, 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 &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 tables; + + /** Headers for tables and output. Contains the name of the solution + components. */ + std::vector headers; + + /** Headers for latex tables. Contains the name of the solution + components. */ + std::vector latex_headers; + + /** Captions for latex. */ + std::vector latex_captions; + + /** Names of the tables. */ + std::vector 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 > data_out; +}; + +#endif diff --git a/deal.II/examples/step-41/include/output_processor.templates.h b/deal.II/examples/step-41/include/output_processor.templates.h new file mode 100644 index 0000000000..06549665a5 --- /dev/null +++ b/deal.II/examples/step-41/include/output_processor.templates.h @@ -0,0 +1,135 @@ +#include "../include/output_processor.h" + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +template +OutputProcessor::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 +void OutputProcessor::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::declare_parameters(prm); + prm.set("Output format", "vtk"); + prm.leave_subsection(); +} +template +void OutputProcessor::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 +void OutputProcessor::prepare_data_output(const DoFHandler &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 partition_int (dh.get_tria().n_active_cells()); + GridTools::get_subdomain_association (dh.get_tria(), partition_int); + Vector partitioning(partition_int.begin(), + partition_int.end()); + static Vector static_partitioning; + static_partitioning.swap(partitioning); + data_out.add_data_vector (static_partitioning, "partitioning"); + } + } + } +} + + +template +void OutputProcessor::add_data_vector(const VECTOR &data_vector, + const std::string &desc) +{ + AssertThrow(initialized, ExcNotInitialized()); + deallog.push("AddingData"); + std::vector 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 +void OutputProcessor::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 +void OutputProcessor::dump_vector (const VECTOR &, + const std::string & ) +{ + Assert(false, ExcNotImplemented()); + // Only specializations exist. +} diff --git a/deal.II/examples/step-41/include/parsed_symmetric_tensor_function.h b/deal.II/examples/step-41/include/parsed_symmetric_tensor_function.h new file mode 100644 index 0000000000..8daeae48b5 --- /dev/null +++ b/deal.II/examples/step-41/include/parsed_symmetric_tensor_function.h @@ -0,0 +1,193 @@ +/* + * Immersed Boundary Problem: + * + * Header Files + * + * Author: + * Luca Heltai + * ============= + * 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 +#include +#include +#include +#include +#include + +#include + +using namespace dealii; + +template +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 & operator() (const Point &p) const; + + /** Get a vector of the values at a number of specified points */ + void value_list(const std::vector< Point > &points, + std::vector< SymmetricTensor > &values); + + /** Get the tensor currently stored. */ + inline const SymmetricTensor & 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 f; + mutable SymmetricTensor t; +}; + + +template +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 +inline +double operator*(const Tensor<2,dim> &gradv, + const SymmetricTensor<2,dim> & C) { + double T = 0; + for(unsigned int i=0; i +inline +double double_contract(const SymmetricTensor<2,dim> & C, + const Tensor<2,dim> &gradv) { + double T = 0; + for(unsigned int i=0; i +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 */ +/* inline */ +/* double double_contract(const Tensor<2,dim> &gradv, */ +/* const Tensor<2,dim> & C) { */ +/* double T = 0; */ +/* for(unsigned int i=0; i */ +/* 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 */ +/* inline */ +/* double operator*(const Tensor<2,dim> &gradv, */ +/* const SymmetricTensor<2,dim> & C) { */ +/* double T = 0; */ +/* for(unsigned int i=0; i */ +/* inline */ +/* double double_contract(const SymmetricTensor<2,dim> & C, */ +/* const Tensor<2,dim> &gradv) { */ +/* double T = 0; */ +/* for(unsigned int i=0; i */ +/* 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 +inline +double double_contract(const Tensor<2,dim> &gradv, + const Tensor<2,dim> & C) { + double T = 0; + for(unsigned int i=0; i + * ============= + * 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 + +template + void ParsedSymmetricTensorFunction::value_list(const std::vector< Point > &points, + std::vector< SymmetricTensor > &values) +{ + + for (unsigned int n=0; noperator()(points[n]); + +} + +template + unsigned int ParsedSymmetricTensorFunction::get_dim_triangular() +{ + if (rank == 2) + return dim; + else if (rank == 4) + return (dim*(dim+1)/2); + else { + AssertThrow(false, ExcInternalError()); + return 0; + } +} + + +template +void ParsedSymmetricTensorFunction::set_time (double t) { + f.set_time(t); +} + + +template +ParsedSymmetricTensorFunction::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 +void ParsedSymmetricTensorFunction::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 expr(get_dim_triangular()); + for(unsigned int i=0; i + void ParsedSymmetricTensorFunction::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 const_list = + Utilities::split_string_list(constants_list, ','); + std::map constants; + for(unsigned int i = 0; i < const_list.size(); ++i) { + std::vector 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.")); + } +} diff --git a/deal.II/examples/step-41/include/utilities.h b/deal.II/examples/step-41/include/utilities.h new file mode 100644 index 0000000000..7ff0f0f6c2 --- /dev/null +++ b/deal.II/examples/step-41/include/utilities.h @@ -0,0 +1,62 @@ +#ifndef UTILITIES_HLT +#define UTILITIES_HLT + +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + +template +class InverseMatrix : public Subscriptor +{ + public: + InverseMatrix (const Matrix &m); + + void vmult (Vector &dst, + const Vector &src) const; + + void Tvmult (Vector &dst, + const Vector &src) const; + + private: + const SmartPointer matrix; + + mutable GrowingVectorMemory<> vector_memory; +}; + + +class BBt : public Subscriptor +{ + public: + void reinit (const BlockSparseMatrix &A); + + + void vmult (Vector &dst, + const Vector &src) const; + + void Tvmult (Vector &dst, + const Vector &src) const; + + private: + SmartPointer > system_matrix; + + mutable BlockVector u, v; + + unsigned int dim; +}; + +template +void smart_delete (SmartPointer &sp) { + if(sp) { + TYPE * p = sp; + sp = 0; + delete p; + } +} + +#endif diff --git a/deal.II/examples/step-41/include/vector_space.h b/deal.II/examples/step-41/include/vector_space.h new file mode 100644 index 0000000000..17b50f9308 --- /dev/null +++ b/deal.II/examples/step-41/include/vector_space.h @@ -0,0 +1,325 @@ +#ifndef VECTOR_SPACE_H +#define VECTOR_SPACE_H + +#include +#include +#include +#include +#include +#include +#include +#include + +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 +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 &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 &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 &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 + 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 target_comps=std::vector()); + + /** 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 + target_components=std::vector()); + + /** 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 & f, + std::map & bvalues); + + /** Return reference to finite element.*/ + inline FiniteElement & get_fe() { + return *fe; + } + + /** Return reference to current dof handler.*/ + inline MGDoFHandler & get_dh() { + return *dh; + } + + /** Return reference to previous dof handler.*/ + inline MGDoFHandler & get_last_dh() { + return *last_dh; + } + + /** Return reference to current triangulation..*/ + inline Triangulation & get_tria() { + return *tria; + } + + /** Return reference to current triangulation..*/ + inline Triangulation & get_last_tria() { + return *last_tria; + } + + /** Return reference to coarse triangulation..*/ + inline Triangulation & get_coarse_tria() { + return *coarse; + } + + /** Return reference to hanging node constraints..*/ + inline ConstraintMatrix & get_hang() { + return hang; + } + + /** Return reference to mapping.*/ + inline Mapping & get_mapping() { + return *mapping; + } + + + /** Return constant reference to finite element.*/ + inline const FiniteElement & get_fe() const { + return *fe; + } + + /** Return constant reference to current dof handler.*/ + inline const MGDoFHandler & get_dh() const { + return *dh; + } + + /** Return constant reference to previous dof handler.*/ + inline const MGDoFHandler & get_last_dh() const { + return *last_dh; + } + + /** Return reference to current triangulation..*/ + inline const Triangulation & get_tria() const { + return *tria; + } + + /** Return reference to coarse triangulation..*/ + inline const Triangulation & get_coarse_tria() const { + return *coarse; + } + + /** Return reference to current triangulation..*/ + inline const Triangulation & 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 & 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 & n_dofs_pp() const { + return local_dofs_per_process; + } + + /** Number of dofs per block. */ + inline const std::vector & 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 > & n_dofs_pb_pp() const { + return dofs_per_block_per_process; + } + + /** Number of dofs per process per block. */ + inline const std::vector > & 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 > dirichlet_bc; + + /** Neumann Boundary Indicators. */ + std::map > neumann_bc; + + /** Other Boundary Indicators. */ + std::map > other_bc; +private: + + /** For internal use. Counts the dofs per block, block and what not. */ + void count_dofs(std::vector target_components); + + std::string fe_name; + std::vector ordering; + unsigned int map_type; + + std::vector d_bc; + std::vector n_bc; + std::vector o_bc; + + /** Pointer to the finite element used. */ + SmartPointer > fe; + + /** Pointer to the dofhandler used */ + SmartPointer > dh; + + /** Pointer to the last dofhandler used */ + SmartPointer > last_dh; + + /** Pointer to a pristine coarse triangulation. */ + SmartPointer > coarse; + + /** Pointer to the current triangulation used */ + SmartPointer > tria; + + /** Pointer to the previous triangulation used. */ + SmartPointer > 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; + + /** 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 local_dofs_per_process; + + /** Local size of vectors. */ + std::vector dofs_per_block; + + /** Local sizes of vectors divided by blocks. */ + std::vector > dofs_per_block_per_process; + + /** Local sizes of vectors divided by blocks. */ + std::vector > 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 wind; +}; +#endif diff --git a/deal.II/examples/step-41/include/vector_space.templates.h b/deal.II/examples/step-41/include/vector_space.templates.h new file mode 100644 index 0000000000..3b62445756 --- /dev/null +++ b/deal.II/examples/step-41/include/vector_space.templates.h @@ -0,0 +1,583 @@ +#include "../include/vector_space.h" +#include "../include/utilities.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +template +VectorSpace::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 +VectorSpace::VectorSpace(ParameterHandler &prm, + Triangulation &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 +VectorSpace::~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 +void VectorSpace::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 +void VectorSpace::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 +void VectorSpace::refine_globally() +{ + deallog.push("GlobalRefinement"); + tria->refine_global (1); + redistribute_dofs(); + deallog.pop(); +} + +template +void VectorSpace::flag_for_refinement(Vector &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 +void VectorSpace::reorder(int level, + std::vector 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&>(*dh)); + } + + for(unsigned int i=0; i&>(*dh)); + } else if(ordering[i] == "comp") { + deallog << "Renumbering component wise." << std::endl; + DoFRenumbering::component_wise(static_cast&>(*dh), target_comps); + } else if(ordering[i] == "upwind") { + DoFRenumbering::downstream(static_cast&>(*dh), wind); + } else if(ordering[i] == "none") {} + else { + AssertThrow(false, ExcMessage("This Reordering not implemented")); + } + } + +} + +template +void VectorSpace::measure_mesh() { + typename Triangulation::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 +void VectorSpace::reinit(ParameterHandler &prm, + Triangulation &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 * fe_p=0; + + try { + fe_p = FETools::get_fe_from_name(fe_name); + } catch (...) { + + // Did not recognize the finite element to use. Try your own + // one... + // fe_p = new FESystem(FE_BGH(2),1, FE_DGP(1), 1); + throw; + } + fe = fe_p; + deallog << "Finite Element Space: " << fe->get_name() << endl; + + if(map_type == 0) { + mapping= new MappingCartesian(); + } else { + mapping = new MappingQ(map_type); + } + + // Store a pointer to the triangulation + coarse = ⅆ + + // Now generate an empty triangulation. + last_tria = new Triangulation(); + + // And a copy of the pristine version of the triangulation. + tria = new Triangulation(); + tria->copy_triangulation(*coarse); + + // And we initialize the current triangulation. + initialize_mesh(); + + // Generate new DoFHandlers + dh = new MGDoFHandler(*tria); + last_dh = new MGDoFHandler(*last_tria); + + // Now generates the boundary maps + std::vector > bcs; + bcs.push_back(d_bc); + bcs.push_back(n_bc); + bcs.push_back(o_bc); + + std::vector > > bc_maps(3); + + for(unsigned int bcnumber=0; bcnumber<3; ++bcnumber) { + // Get some aliases for the bcs + std::vector & bc = bcs[bcnumber]; + std::map > & bc_map = bc_maps[bcnumber]; + + if(bc.size()) + for(unsigned int i=0; i < bc.size(); ++i) { + std::vector id_and_comps = + Utilities::split_string_list(bc[i], ':'); + AssertThrow(id_and_comps.size() == 2, + ExcMessage("Wrong Format for boundary indicator map.")); + std::vector ids = + Utilities::string_to_int(Utilities::split_string_list(id_and_comps[0])); + std::vector comps = + Utilities::string_to_int(Utilities::split_string_list(id_and_comps[1])); + + unsigned int n_c = fe->n_components(); + std::vector filter(n_c, false); + // Now check that the components make sense + for(unsigned int i=0; i +void VectorSpace::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 +void VectorSpace::initialize_mesh() +{ + deallog.push("InitializeMesh"); + tria->refine_global (global_refinement); + deallog << "Active Cells: " + << tria->n_active_cells() + << endl; + deallog.pop(); +} + +template +void VectorSpace::redistribute_dofs(std::vector 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 +void VectorSpace::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 +void VectorSpace::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 wind_str = + Utilities::split_string_list(prm.get ("Wind direction") ); + for(unsigned int i=0; (i +void VectorSpace::interpolate_dirichlet_bc(const Function & f, + std::map & bvalues) +{ + deallog.push("DBC"); + std::map >::iterator + dmap = dirichlet_bc.begin(), + dmapend = dirichlet_bc.end(); + + unsigned int last_counted = 0; + for(; dmap!=dmapend; ++dmap) { + char id = dmap->first; + std::vector &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&) *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 +void VectorSpace::count_dofs(std::vector 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 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(n_mpi_processes, 0)); + dofs_per_process_per_block.resize(n_mpi_processes, + std::vector(number_of_blocks, 0)); + + // We now have a vector that contains all subdomain ids for each dof + std::vector 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 mask(fe->n_components(), false); + for(unsigned int i=0; i this_block_dofs(dh->n_dofs()); + DoFTools::extract_dofs(*dh, mask, this_block_dofs); + + // Reset counter + for(unsigned int j=0; jn_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; jn_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(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]; +} diff --git a/deal.II/examples/step-41/mesh/square.inp b/deal.II/examples/step-41/mesh/square.inp new file mode 100644 index 0000000000..23f812759a --- /dev/null +++ b/deal.II/examples/step-41/mesh/square.inp @@ -0,0 +1,16 @@ +# 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 diff --git a/deal.II/examples/step-41/source/base.cc b/deal.II/examples/step-41/source/base.cc new file mode 100644 index 0000000000..9da975a53b --- /dev/null +++ b/deal.II/examples/step-41/source/base.cc @@ -0,0 +1,3 @@ +#include "../include/base.templates.h" + +template class Base; diff --git a/deal.II/examples/step-41/source/camclay.cc b/deal.II/examples/step-41/source/camclay.cc new file mode 100644 index 0000000000..506ae0cee0 --- /dev/null +++ b/deal.II/examples/step-41/source/camclay.cc @@ -0,0 +1,3 @@ +#include "../include/camclay.templates.h" + +template class CamClay; diff --git a/deal.II/examples/step-41/source/domain.cc b/deal.II/examples/step-41/source/domain.cc new file mode 100644 index 0000000000..fa3bc745c3 --- /dev/null +++ b/deal.II/examples/step-41/source/domain.cc @@ -0,0 +1,3 @@ +#include "../include/domain.templates.h" + +template class Domain; diff --git a/deal.II/examples/step-41/source/error_handler.cc b/deal.II/examples/step-41/source/error_handler.cc new file mode 100644 index 0000000000..0403a7cfd2 --- /dev/null +++ b/deal.II/examples/step-41/source/error_handler.cc @@ -0,0 +1,11 @@ +#include "../include/error_handler.templates.h" +#include +#include +#include + +template class ErrorHandler >; +template class ErrorHandler >; +#ifdef CONTRIB_USE_PETSC +template class ErrorHandler; +#endif +//template class ErrorHandler; diff --git a/deal.II/examples/step-41/source/linear_elastic.cc b/deal.II/examples/step-41/source/linear_elastic.cc new file mode 100644 index 0000000000..f90ad9b616 --- /dev/null +++ b/deal.II/examples/step-41/source/linear_elastic.cc @@ -0,0 +1,4 @@ +#include "../include/linear_elastic.templates.h" + + +template class LinearElastic; diff --git a/deal.II/examples/step-41/source/local_assemble.cc b/deal.II/examples/step-41/source/local_assemble.cc new file mode 100644 index 0000000000..54ac66cf2a --- /dev/null +++ b/deal.II/examples/step-41/source/local_assemble.cc @@ -0,0 +1,6 @@ +#include "../include/local_assemble.templates.h" +#include +#include + +template class LocalAssemble >; +template class LocalAssemble >; diff --git a/deal.II/examples/step-41/source/local_assemble_base.cc b/deal.II/examples/step-41/source/local_assemble_base.cc new file mode 100644 index 0000000000..5ab1fc3c62 --- /dev/null +++ b/deal.II/examples/step-41/source/local_assemble_base.cc @@ -0,0 +1,5 @@ +#include "../include/local_assemble_base.templates.h" + +template class LocalAssembleBase >; +template class LocalAssembleBase >; + diff --git a/deal.II/examples/step-41/source/local_assemble_elastic_matrix.cc b/deal.II/examples/step-41/source/local_assemble_elastic_matrix.cc new file mode 100644 index 0000000000..1b704c5f51 --- /dev/null +++ b/deal.II/examples/step-41/source/local_assemble_elastic_matrix.cc @@ -0,0 +1,4 @@ +#include "../include/local_assemble_elastic_matrix.templates.h" + + +template class LocalAssembleElasticMatrix; diff --git a/deal.II/examples/step-41/source/local_assemble_elastic_rhs.cc b/deal.II/examples/step-41/source/local_assemble_elastic_rhs.cc new file mode 100644 index 0000000000..0c351d4be1 --- /dev/null +++ b/deal.II/examples/step-41/source/local_assemble_elastic_rhs.cc @@ -0,0 +1,4 @@ +#include "../include/local_assemble_elastic_rhs.templates.h" + + +template class LocalAssembleElasticRHS; diff --git a/deal.II/examples/step-41/source/local_assemble_mass.cc b/deal.II/examples/step-41/source/local_assemble_mass.cc new file mode 100644 index 0000000000..78d0dc9954 --- /dev/null +++ b/deal.II/examples/step-41/source/local_assemble_mass.cc @@ -0,0 +1,6 @@ +#include "../include/local_assemble_mass.templates.h" +#include +#include + +template class LocalAssembleMass >; +template class LocalAssembleMass >; diff --git a/deal.II/examples/step-41/source/local_assemble_plastic_project.cc b/deal.II/examples/step-41/source/local_assemble_plastic_project.cc new file mode 100644 index 0000000000..a7600a71f9 --- /dev/null +++ b/deal.II/examples/step-41/source/local_assemble_plastic_project.cc @@ -0,0 +1,4 @@ +#include "../include/local_assemble_plastic_project.templates.h" + + +template class LocalAssemblePlasticProject; diff --git a/deal.II/examples/step-41/source/local_assemble_scalar_project.cc b/deal.II/examples/step-41/source/local_assemble_scalar_project.cc new file mode 100644 index 0000000000..c8075e678d --- /dev/null +++ b/deal.II/examples/step-41/source/local_assemble_scalar_project.cc @@ -0,0 +1,4 @@ +#include "../include/local_assemble_scalar_project.templates.h" + + +template class LocalAssembleScalarProject; diff --git a/deal.II/examples/step-41/source/main.cc b/deal.II/examples/step-41/source/main.cc new file mode 100644 index 0000000000..e45e1b592f --- /dev/null +++ b/deal.II/examples/step-41/source/main.cc @@ -0,0 +1,37 @@ +#include +#include "../include/base.h" + +int main () +{ + try + { + Base 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; +} diff --git a/deal.II/examples/step-41/source/output_processor.cc b/deal.II/examples/step-41/source/output_processor.cc new file mode 100644 index 0000000000..dd5ff0027e --- /dev/null +++ b/deal.II/examples/step-41/source/output_processor.cc @@ -0,0 +1,26 @@ +#include "../include/output_processor.templates.h" +#include +#include +#include + +template <> +void OutputProcessor >::dump_vector (const Vector &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 >; + +template class OutputProcessor >; +template class OutputProcessor >; + +#ifdef DEAL_II_USE_PETSC +template class OutputProcessor; +#endif + +// template class OutputProcessor; diff --git a/deal.II/examples/step-41/source/parsed_symmetric_tensor_function.cc b/deal.II/examples/step-41/source/parsed_symmetric_tensor_function.cc new file mode 100644 index 0000000000..3c1df0507c --- /dev/null +++ b/deal.II/examples/step-41/source/parsed_symmetric_tensor_function.cc @@ -0,0 +1,158 @@ +/* + * Immersed Boundary Problem: + * + * Header Files + * + * Author: + * Luca Heltai + * ============= + * 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>; diff --git a/deal.II/examples/step-41/source/vector_space.cc b/deal.II/examples/step-41/source/vector_space.cc new file mode 100644 index 0000000000..fa06236c7f --- /dev/null +++ b/deal.II/examples/step-41/source/vector_space.cc @@ -0,0 +1,3 @@ +#include "../include/vector_space.templates.h" + +template class VectorSpace; -- 2.39.5