From: heltai Date: Tue, 11 Oct 2011 11:53:12 +0000 (+0000) Subject: Removed unfinished step-41 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bd679d5ff9449bfd33c8a7e7e5f291b871b62afc;p=dealii-svn.git Removed unfinished step-41 git-svn-id: https://svn.dealii.org/trunk@24591 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-1/Makefile b/deal.II/examples/step-1/Makefile index 084f822e9a..b4c21ed325 100644 --- a/deal.II/examples/step-1/Makefile +++ b/deal.II/examples/step-1/Makefile @@ -33,8 +33,6 @@ D = ../../ clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2 - - # # # Usually, you will not need to change anything beyond this point. @@ -45,6 +43,41 @@ clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2 # settings include $D/common/Make.global_options +CXX += -I$(VTK_DIR)/include/vtk-5.6 +LIBS += -L/usr/local/vtk/lib/vtk/ -I$(VTK_DIR)/include/vtk-5.6 \ +-lvtkCharts \ +-lvtkCommon \ +-lvtkDICOMParser \ +-lvtkFiltering \ +-lvtkGenericFiltering \ +-lvtkGeovis \ +-lvtkGraphics \ +-lvtkHybrid \ +-lvtkIO \ +-lvtkImaging \ +-lvtkInfovis \ +-lvtkNetCDF \ +-lvtkRendering \ +-lvtkViews \ +-lvtkVolumeRendering \ +-lvtkWidgets \ +-lvtkalglib \ +-lvtkexoIIc \ +-lvtkexpat \ +-lvtkfreetype \ +-lvtkftgl \ +-lvtkjpeg \ +-lvtklibxml2 \ +-lvtkmetaio \ +-lvtkpng \ +-lvtkproj4 \ +-lvtksqlite \ +-lvtksys \ +-lvtktiff \ +-lvtkverdict \ +-lvtkzlib + + # Since the whole project consists of only one file, we need not # consider difficult dependencies. We only have to declare the diff --git a/deal.II/examples/step-34/parameters.prm b/deal.II/examples/step-34/parameters.prm index 5e12caf422..03e374ec0a 100644 --- a/deal.II/examples/step-34/parameters.prm +++ b/deal.II/examples/step-34/parameters.prm @@ -1,10 +1,10 @@ # Listing of Parameters # --------------------- set Extend solution on the -2,2 box = true -set External refinement = 5 +set External refinement = 7 set Number of cycles = 4 set Run 2d simulation = true -set Run 3d simulation = true +set Run 3d simulation = false subsection Exact solution 2d diff --git a/deal.II/examples/step-41/CamClay.prm b/deal.II/examples/step-41/CamClay.prm deleted file mode 100644 index 404fbd435c..0000000000 --- a/deal.II/examples/step-41/CamClay.prm +++ /dev/null @@ -1,306 +0,0 @@ -# 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 deleted file mode 100644 index 6ea309f791..0000000000 --- a/deal.II/examples/step-41/Makefile +++ /dev/null @@ -1,149 +0,0 @@ -# $Id$ - -# The large projects Makefile looks much like the one for small -# projects. Basically, only the following seven parameters need to be -# set by you: - -application-name = camclay - -# The first denotes the dimension for which the program is to be -# compiled: -deal_II_dimension = 2 - -# The second tells us the name of the executable. It is prefixed by -# `lib/' to designate its destination directory. Note that the program -# name depends on the dimension, so you can keep copies for the -# different dimensions around: -target = lib/$(application-name)-$(deal_II_dimension)d - -# The `debug-mode' variable works as in the small projects Makefile: -debug-mode = on - -# And so does the following variable. You will have to set it to -# something more reasonable, of course. -clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2 - -# Finally, here is a variable which tells the `run' rule which -# parameters to pass to the executable. Usually, this will be the name -# of an input file. -run-parameters = parameter-file.prm - -# Now, this is the last variable you need to set, namely the path to -# the deal.II toplevel directory: -D = ../../ - -# -# -# Usually, you will not need to change anything beyond this point. -# -# -# This tells `make' where to find the global settings and rules: -include $D/common/Make.global_options - - -################################################################ -# This example program will only work if Trilinos is installed. If this -# is not the case, then simply redefine the main targets to do nothing -ifneq ($(USE_CONTRIB_TRILINOS),yes) -default run clean: - @echo - @echo "===========================================================" - @echo "= This program cannot be compiled without Trilinos. Make=" - @echo "= sure you have Trilinos installed and detected during =" - @echo "= configuration of deal.II =" - @echo "===========================================================" - @echo - -else -# -################################################################ - - -# First get a list of files belonging to the project. Include files -# are expected in `include/', while implementation files are expected -# in `source/'. Object files are placed into `lib/[123]d', using the -# same base name as the `.cc' file. -cc-files = $(shell echo source/*.cc) -o-files = $(cc-files:source/%.cc=lib/$(deal_II_dimension)d/%.$(OBJEXT)) -go-files = $(cc-files:source/%.cc=lib/$(deal_II_dimension)d/%.g.$(OBJEXT)) -h-files = $(wildcard include/*.h) -lib-h-files = $(shell echo $D/include/deal.II/*/*.h) - -# Since the whole project consists of only one file, we need not -# consider difficult dependencies. We only have to declare the -# libraries which we want to link to the object file. deal.II has two -# libraries: one for the debug mode version of the -# application and one for optimized mode. -libs.g := $(lib-deal2.g) -libs.o := $(lib-deal2.o) - - - -# Now use the information from above to define the set of libraries to -# link with and the flags to be passed to the compiler: -ifeq ($(debug-mode),on) - libraries = $(go-files) $(libs.g) - flags = $(CXXFLAGS.g) -else - libraries = $(o-files) $(libs.o) - flags = $(CXXFLAGS.o) -endif - - -# Then augment the compiler flags by a specification of the dimension -# for which the program shall be compiled: -flags += -Ddeal_II_dimension=$(deal_II_dimension) - - -# The following two rules define how to compile C++ files into object -# files: -lib/$(deal_II_dimension)d/%.g.$(OBJEXT) : - @echo =====$(application-name)=======$(deal_II_dimension)d====debug=====$(MT)== $( $@ - -include lib/Makefile.dep - -endif # USE_CONTRIB_TRILINOS diff --git a/deal.II/examples/step-41/doc/builds-on b/deal.II/examples/step-41/doc/builds-on deleted file mode 100644 index 8b13789179..0000000000 --- a/deal.II/examples/step-41/doc/builds-on +++ /dev/null @@ -1 +0,0 @@ - diff --git a/deal.II/examples/step-41/doc/intro.dox b/deal.II/examples/step-41/doc/intro.dox deleted file mode 100644 index 446bf8b910..0000000000 --- a/deal.II/examples/step-41/doc/intro.dox +++ /dev/null @@ -1,4 +0,0 @@ - -

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 deleted file mode 100644 index 56e049c91a..0000000000 --- a/deal.II/examples/step-41/doc/kind +++ /dev/null @@ -1 +0,0 @@ -solids diff --git a/deal.II/examples/step-41/doc/results.dox b/deal.II/examples/step-41/doc/results.dox deleted file mode 100644 index 6121fb67ed..0000000000 --- a/deal.II/examples/step-41/doc/results.dox +++ /dev/null @@ -1,3 +0,0 @@ -

Results

- -Astonishing results here. diff --git a/deal.II/examples/step-41/doc/tooltip b/deal.II/examples/step-41/doc/tooltip deleted file mode 100644 index 03e9264e3d..0000000000 --- a/deal.II/examples/step-41/doc/tooltip +++ /dev/null @@ -1 +0,0 @@ -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 deleted file mode 100644 index 544f4322e6..0000000000 --- a/deal.II/examples/step-41/include/assemble_flags.h +++ /dev/null @@ -1,122 +0,0 @@ -//--------------------------------------------------------------------------- -// $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 deleted file mode 100644 index 1a13228386..0000000000 --- a/deal.II/examples/step-41/include/base.h +++ /dev/null @@ -1,113 +0,0 @@ -#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 deleted file mode 100644 index 5f6f4c63ed..0000000000 --- a/deal.II/examples/step-41/include/base.templates.h +++ /dev/null @@ -1,468 +0,0 @@ -#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.size(), - std::vector >(dim)); - std::vector< Point > points(qf.size()); - - 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.size(), - std::vector >(dim)); - std::vector< std::vector< Tensor<1,dim> > > plastic_grads(qf.size(), - std::vector >(dim)); - std::vector< Vector > hard_values(qf.size(), Vector(dim)); - std::vector< Point > points(qf.size()); - - //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 deleted file mode 100644 index e1ec25367b..0000000000 --- a/deal.II/examples/step-41/include/camclay.h +++ /dev/null @@ -1,172 +0,0 @@ -#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 deleted file mode 100644 index ca9c2bad9e..0000000000 --- a/deal.II/examples/step-41/include/camclay.templates.h +++ /dev/null @@ -1,670 +0,0 @@ -#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.size(); - - - 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.size(); - - 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 deleted file mode 100644 index e4967c0239..0000000000 --- a/deal.II/examples/step-41/include/domain.h +++ /dev/null @@ -1,72 +0,0 @@ -#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 deleted file mode 100644 index 188f006e48..0000000000 --- a/deal.II/examples/step-41/include/domain.templates.h +++ /dev/null @@ -1,120 +0,0 @@ -#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 deleted file mode 100644 index fbc52b6c23..0000000000 --- a/deal.II/examples/step-41/include/error_handler.h +++ /dev/null @@ -1,177 +0,0 @@ -#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 deleted file mode 100644 index 2b646b649b..0000000000 --- a/deal.II/examples/step-41/include/error_handler.templates.h +++ /dev/null @@ -1,396 +0,0 @@ -#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 deleted file mode 100644 index 02b51cf708..0000000000 --- a/deal.II/examples/step-41/include/linear_elastic.templates.h +++ /dev/null @@ -1,160 +0,0 @@ -#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 deleted file mode 100644 index 89c0bdb215..0000000000 --- a/deal.II/examples/step-41/include/local_assemble.h +++ /dev/null @@ -1,66 +0,0 @@ -#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 deleted file mode 100644 index 1a2d4a598e..0000000000 --- a/deal.II/examples/step-41/include/local_assemble.templates.h +++ /dev/null @@ -1,67 +0,0 @@ -#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 deleted file mode 100644 index d8576f33ef..0000000000 --- a/deal.II/examples/step-41/include/local_assemble_base.h +++ /dev/null @@ -1,102 +0,0 @@ -#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 deleted file mode 100644 index 98fdc4e873..0000000000 --- a/deal.II/examples/step-41/include/local_assemble_base.templates.h +++ /dev/null @@ -1,68 +0,0 @@ -#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 deleted file mode 100644 index 64c21f2855..0000000000 --- a/deal.II/examples/step-41/include/local_assemble_elastic_matrix.h +++ /dev/null @@ -1,50 +0,0 @@ -#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 deleted file mode 100644 index 76480c4375..0000000000 --- a/deal.II/examples/step-41/include/local_assemble_elastic_matrix.templates.h +++ /dev/null @@ -1,85 +0,0 @@ -#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 deleted file mode 100644 index a31081ee43..0000000000 --- a/deal.II/examples/step-41/include/local_assemble_elastic_rhs.h +++ /dev/null @@ -1,65 +0,0 @@ -#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 deleted file mode 100644 index c428b2671b..0000000000 --- a/deal.II/examples/step-41/include/local_assemble_elastic_rhs.templates.h +++ /dev/null @@ -1,133 +0,0 @@ -#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 deleted file mode 100644 index fc4be24031..0000000000 --- a/deal.II/examples/step-41/include/local_assemble_hardening_project.h +++ /dev/null @@ -1,53 +0,0 @@ -#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 deleted file mode 100644 index 9d3e8085cc..0000000000 --- a/deal.II/examples/step-41/include/local_assemble_mass.h +++ /dev/null @@ -1,45 +0,0 @@ -#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 deleted file mode 100644 index da022244eb..0000000000 --- a/deal.II/examples/step-41/include/local_assemble_mass.templates.h +++ /dev/null @@ -1,93 +0,0 @@ -#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 deleted file mode 100644 index 5bd80c81f6..0000000000 --- a/deal.II/examples/step-41/include/local_assemble_plastic_project.h +++ /dev/null @@ -1,53 +0,0 @@ -#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 deleted file mode 100644 index 3105ddeb58..0000000000 --- a/deal.II/examples/step-41/include/local_assemble_plastic_project.templates.h +++ /dev/null @@ -1,93 +0,0 @@ -#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 deleted file mode 100644 index 6c05742d2d..0000000000 --- a/deal.II/examples/step-41/include/local_assemble_scalar_project.h +++ /dev/null @@ -1,51 +0,0 @@ -#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 deleted file mode 100644 index 3821b90480..0000000000 --- a/deal.II/examples/step-41/include/local_assemble_scalar_project.templates.h +++ /dev/null @@ -1,84 +0,0 @@ -#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 deleted file mode 100644 index dd80259187..0000000000 --- a/deal.II/examples/step-41/include/my_tools.h +++ /dev/null @@ -1,989 +0,0 @@ -//--------------------------------------------------------------------------- -// $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 deleted file mode 100644 index daabff1d70..0000000000 --- a/deal.II/examples/step-41/include/output_processor.templates.h +++ /dev/null @@ -1,135 +0,0 @@ -#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 deleted file mode 100644 index 6a93c0bd7e..0000000000 --- a/deal.II/examples/step-41/include/parsed_symmetric_tensor_function.h +++ /dev/null @@ -1,193 +0,0 @@ -/* - * 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 deleted file mode 100644 index e1cdf835fd..0000000000 --- a/deal.II/examples/step-41/include/utilities.h +++ /dev/null @@ -1,62 +0,0 @@ -#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 deleted file mode 100644 index 1d200f7abd..0000000000 --- a/deal.II/examples/step-41/include/vector_space.h +++ /dev/null @@ -1,325 +0,0 @@ -#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 deleted file mode 100644 index 1a10fb3adf..0000000000 --- a/deal.II/examples/step-41/include/vector_space.templates.h +++ /dev/null @@ -1,583 +0,0 @@ -#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 deleted file mode 100644 index 23f812759a..0000000000 --- a/deal.II/examples/step-41/mesh/square.inp +++ /dev/null @@ -1,16 +0,0 @@ -# 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 deleted file mode 100644 index 9da975a53b..0000000000 --- a/deal.II/examples/step-41/source/base.cc +++ /dev/null @@ -1,3 +0,0 @@ -#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 deleted file mode 100644 index 506ae0cee0..0000000000 --- a/deal.II/examples/step-41/source/camclay.cc +++ /dev/null @@ -1,3 +0,0 @@ -#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 deleted file mode 100644 index fa3bc745c3..0000000000 --- a/deal.II/examples/step-41/source/domain.cc +++ /dev/null @@ -1,3 +0,0 @@ -#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 deleted file mode 100644 index f4a7e5ad58..0000000000 --- a/deal.II/examples/step-41/source/error_handler.cc +++ /dev/null @@ -1,11 +0,0 @@ -#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 deleted file mode 100644 index f90ad9b616..0000000000 --- a/deal.II/examples/step-41/source/linear_elastic.cc +++ /dev/null @@ -1,4 +0,0 @@ -#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 deleted file mode 100644 index 338332869f..0000000000 --- a/deal.II/examples/step-41/source/local_assemble.cc +++ /dev/null @@ -1,6 +0,0 @@ -#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 deleted file mode 100644 index 5ab1fc3c62..0000000000 --- a/deal.II/examples/step-41/source/local_assemble_base.cc +++ /dev/null @@ -1,5 +0,0 @@ -#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 deleted file mode 100644 index 1b704c5f51..0000000000 --- a/deal.II/examples/step-41/source/local_assemble_elastic_matrix.cc +++ /dev/null @@ -1,4 +0,0 @@ -#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 deleted file mode 100644 index 0c351d4be1..0000000000 --- a/deal.II/examples/step-41/source/local_assemble_elastic_rhs.cc +++ /dev/null @@ -1,4 +0,0 @@ -#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 deleted file mode 100644 index f6d6e9f5ac..0000000000 --- a/deal.II/examples/step-41/source/local_assemble_mass.cc +++ /dev/null @@ -1,6 +0,0 @@ -#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 deleted file mode 100644 index a7600a71f9..0000000000 --- a/deal.II/examples/step-41/source/local_assemble_plastic_project.cc +++ /dev/null @@ -1,4 +0,0 @@ -#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 deleted file mode 100644 index c8075e678d..0000000000 --- a/deal.II/examples/step-41/source/local_assemble_scalar_project.cc +++ /dev/null @@ -1,4 +0,0 @@ -#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 deleted file mode 100644 index e45e1b592f..0000000000 --- a/deal.II/examples/step-41/source/main.cc +++ /dev/null @@ -1,37 +0,0 @@ -#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 deleted file mode 100644 index 2b28d33a9e..0000000000 --- a/deal.II/examples/step-41/source/output_processor.cc +++ /dev/null @@ -1,26 +0,0 @@ -#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 deleted file mode 100644 index 3c1df0507c..0000000000 --- a/deal.II/examples/step-41/source/parsed_symmetric_tensor_function.cc +++ /dev/null @@ -1,158 +0,0 @@ -/* - * 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 deleted file mode 100644 index fa06236c7f..0000000000 --- a/deal.II/examples/step-41/source/vector_space.cc +++ /dev/null @@ -1,3 +0,0 @@ -#include "../include/vector_space.templates.h" - -template class VectorSpace; diff --git a/deal.II/examples/step-41/step-41.cc b/deal.II/examples/step-41/step-41.cc deleted file mode 100644 index e45e1b592f..0000000000 --- a/deal.II/examples/step-41/step-41.cc +++ /dev/null @@ -1,37 +0,0 @@ -#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-6/step-6.cc b/deal.II/examples/step-6/step-6.cc index 657a12ade4..5401616f01 100644 --- a/deal.II/examples/step-6/step-6.cc +++ b/deal.II/examples/step-6/step-6.cc @@ -16,28 +16,30 @@ // been covered in previous examples // and will thus not be further // commented on. -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "constrained_matrix.h" #include #include @@ -51,7 +53,7 @@ // tri-linear elements, but we will // now use it for bi-quadratic // elements: -#include +#include // We will not read the grid from a // file as in the previous example, // but generate it using a function @@ -61,7 +63,7 @@ // the solution) in each step, so we // need the following include file // instead of grid_in.h: -#include +#include // When using locally refined grids, @@ -77,7 +79,7 @@ // following file contains a class // which is used to handle these // constraints: -#include +#include // In order to refine our grids // locally, we need a function from @@ -86,7 +88,7 @@ // coarsening based on the error // indicators we have computed. This // function is defined here: -#include +#include // Finally, we need a simple way to // actually compute the refinement @@ -98,7 +100,7 @@ // often yields quite nicely adapted // grids for a wide class of // problems. -#include +#include // Finally, this is as in previous // programs: @@ -615,8 +617,9 @@ void Step6::assemble_system () // constraints due to hanging // nodes. This is done using the // following two function calls: - hanging_node_constraints.condense (system_matrix); - hanging_node_constraints.condense (system_rhs); + // hanging_node_constraints.condense (system_matrix); + // hanging_node_constraints.condense (system_rhs); + // Using them, degrees of freedom // associated to hanging nodes have // been removed from the linear @@ -695,16 +698,28 @@ void Step6::assemble_system () template void Step6::solve () { - SolverControl solver_control (1000, 1e-12); - SolverCG<> solver (solver_control); + SolverControl solver_control (10000, 1e-12); + SolverGMRES<> solver (solver_control); - PreconditionSSOR<> preconditioner; - preconditioner.initialize(system_matrix, 1.2); + PreconditionIdentity preconditioner; +//preconditioner.initialize(system_matrix, 1.2);rGM - solver.solve (system_matrix, solution, system_rhs, - preconditioner); + double scaling = 10000; + typename DoFHandler::active_cell_iterator + cell=dof_handler.begin_active(), + endc = dof_handler.end(); + for(; cell != endc; ++cell) + scaling = std::min(scaling, cell->diameter()*cell->diameter()); + + ConstrainedOperator, SparseMatrix > + cc(system_matrix, hanging_node_constraints, scaling); - hanging_node_constraints.distribute (solution); + // cc.distribute_rhs(system_rhs); + solution = system_rhs; + cc.distribute_rhs(system_rhs); + + solver.solve (cc, solution, system_rhs, + preconditioner); } @@ -978,6 +993,19 @@ void Step6::output_results (const unsigned int cycle) const GridOut grid_out; grid_out.write_eps (triangulation, output); + + DataOut data_out; + + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.build_patches (); + + filename = "solution-"; + filename += ('0' + cycle); + filename += ".vtu"; + output.close(); + output.open(filename.c_str()); + data_out.write_vtu (output); } @@ -1090,8 +1118,8 @@ void Step6::run () data_out.add_data_vector (solution, "solution"); data_out.build_patches (); - std::ofstream output ("final-solution.eps"); - data_out.write_eps (output); + std::ofstream output ("final-solution.vtk"); + data_out.write_vtk (output); }