From: kanschat Date: Tue, 6 Sep 2011 16:15:40 +0000 (+0000) Subject: start splitting assembler.h X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b42b08f5f5e2f25196079975a585e6dfb804a262;p=dealii-svn.git start splitting assembler.h git-svn-id: https://svn.dealii.org/trunk@24273 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/meshworker/assembler.h b/deal.II/include/deal.II/meshworker/assembler.h index 2f968c3f7e..54acf9779f 100644 --- a/deal.II/include/deal.II/meshworker/assembler.h +++ b/deal.II/include/deal.II/meshworker/assembler.h @@ -18,6 +18,8 @@ #include #include #include +#include +#include #include @@ -78,229 +80,6 @@ namespace MeshWorker */ namespace Assembler { -/** - * The class assembling local contributions to a functional into the - * global functionals. - * - * - * - * @ingroup MeshWorker - * @author Guido Kanschat, 2009 - */ - template - class Functional - { - public: - /** - * Initialize local data to - * store functionals. The - * number n is the - * number of functionals to - * be computed. - */ - void initialize(unsigned int n); - /** - * Initialize the local data - * in the DoFInfo object used - * later for assembling. - * - * The info object refers to - * a cell if - * !face, or - * else to an interior or - * boundary face. - */ - template - void initialize_info(DOFINFO& info, bool face); - - /** - * Assemble the local values - * into the global vectors. - */ - template - void assemble(const DOFINFO& info); - - /** - * Assemble both local values - * into the global vectors. - */ - template - void assemble(const DOFINFO& info1, - const DOFINFO& info2); - - /** - * The value of the ith entry - * in #results. - */ - number operator() (unsigned int i) const; - private: - /** - * The values into which the - * results are added. - */ - std::vector results; - }; - -/** - * Compute cell and face contributions of one or several functionals, - * typically for error estimates. - * - * @ingroup MeshWorker - * @author Guido Kanschat, 2009 - */ - template - class CellsAndFaces - { - public: - /** - * The data type for - * communicating the cell and - * face vectors. - */ - typedef NamedData*> DataVectors; - - /** - * The initialization - * function, specifying the - * #results vectors and - * whether face data should - * be collected separately. - * - * #results should contain - * two block vectors named - * "cells" and "faces" (the - * latter only if - * #separate_faces is - * true). In each of the two, - * each block should have - * equal size and be large - * enough to accomodate all - * user indices set in the - * cells and faces covered by - * the loop it is used - * in. Typically, for - * estimators, this is - * Triangulation::n_active_cells() - * and - * Triangulation::n_faces(), - * respectively. - * - * The use of BlockVector may - * seem cumbersome, but it - * allows us to assemble - * several functionals at the - * same time, one in each - * block. The typical - * situation for error - * estimate is just having a - * single block in each vector. - */ - void initialize(DataVectors& results, - bool separate_faces = true); - /** - * Initialize the local data - * in the - * DoFInfo - * object used later for - * assembling. - * - * The info object refers to - * a cell if - * !face, or - * else to an interior or - * boundary face. - */ - template - void initialize_info(DOFINFO& info, bool face) const; - - /** - * Assemble the local values - * into the global vectors. - */ - template - void assemble(const DOFINFO& info); - - /** - * Assemble both local values - * into the global vectors. - */ - template - void assemble(const DOFINFO& info1, - const DOFINFO& info2); - - /** - * The value of the ith entry - * in #results. - */ - number operator() (unsigned int i) const; - private: - DataVectors results; - bool separate_faces; - }; - - -/** - * Assemble residuals without block structure. - * - * The data structure for this Assembler class is a simple vector on - * each cell with entries from zero to - * FiniteElementData::dofs_per_cell and a simple global vector with - * entries numbered from zero to DoFHandler::n_dofs(). No BlockInfo is - * required and the global vector may be any type of vector having - * element access through operator() (unsigned int) - * - * @ingroup MeshWorker - * @author Guido Kanschat, 2009 - */ - template - class ResidualSimple - { - public: - void initialize(NamedData& results); - /** - * Initialize the constraints. - */ - void initialize(const ConstraintMatrix& constraints); - /** - * Initialize the local data - * in the DoFInfo object used - * later for assembling. - * - * The info object refers to - * a cell if - * !face, or - * else to an interior or - * boundary face. - */ - template - void initialize_info(DOFINFO& info, bool face) const; - - /** - * Assemble the local residuals - * into the global residuals. - */ - template - void assemble(const DOFINFO& info); - - /** - * Assemble both local residuals - * into the global residuals. - */ - template - void assemble(const DOFINFO& info1, - const DOFINFO& info2); - private: - /** - * The global residal vectors - * filled by assemble(). - */ - NamedData > > residuals; - /** - * A pointer to the object containing constraints. - */ - SmartPointer > constraints; - }; - /** * Assemble local residuals into global residuals. * @@ -343,373 +122,61 @@ namespace MeshWorker * DoFInfo * object used later for * assembling. - * - * The info object refers to - * a cell if - * !face, or - * else to an interior or - * boundary face. - */ - template - void initialize_info(DOFINFO& info, bool face) const; - - - /** - * Assemble the local residuals - * into the global residuals. - */ - template - void assemble(const DOFINFO& info); - - /** - * Assemble both local residuals - * into the global residuals. - */ - template - void assemble(const DOFINFO& info1, - const DOFINFO& info2); - private: - /** - * Assemble a single local - * residual into the global. - */ - void assemble(VECTOR& global, - const BlockVector& local, - const std::vector& dof); - - /** - * The global matrices, - * stored as a vector of - * pointers. - */ - NamedData > > residuals; - - /** - * A pointer to the object containing the block structure. - */ - SmartPointer > block_info; - /** - * A pointer to the object containing constraints. - */ - SmartPointer > constraints; - }; - - -/** - * Assemble local matrices into a single global matrix without using - * block structure. - * - * After being initialized with a SparseMatrix object (or another - * matrix offering the same functionality as SparseMatrix::add()), - * this class can be used in a MeshWorker::loop() to assemble the cell - * and face matrices into the global matrix. - * - * @todo On locally refined meshes, a ConstraintMatrix should be used - * to automatically eliminate hanging nodes. - * - * @ingroup MeshWorker - * @author Guido Kanschat, 2009 - */ - template - class MatrixSimple - { - public: - /** - * Constructor, initializing - * the #threshold, which - * limits how small numbers - * may be to be entered into - * the matrix. - */ - MatrixSimple(double threshold = 1.e-12); - - /** - * Store the result matrix - * for later assembling. - */ - void initialize(MATRIX& m); - /** - * Initialize the constraints. - */ - void initialize(const ConstraintMatrix& constraints); - /** - * Initialize the local data - * in the - * DoFInfo - * object used later for - * assembling. - * - * The info object refers to - * a cell if - * !face, or - * else to an interior or - * boundary face. - */ - template - void initialize_info(DOFINFO& info, bool face) const; - - /** - * Assemble the matrix - * DoFInfo::M1[0] - * into the global matrix. - */ - template - void assemble(const DOFINFO& info); - - /** - * Assemble both local - * matrices in the info - * objects into the global - * matrix. - */ - template - void assemble(const DOFINFO& info1, - const DOFINFO& info2); - private: - /** - * Assemble a single matrix - * into #matrix. - */ - void assemble(const FullMatrix& M, - const std::vector& i1, - const std::vector& i2); - - /** - * The global matrix being - * assembled. - */ - SmartPointer > matrix; - /** - * A pointer to the object containing constraints. - */ - SmartPointer > constraints; - - /** - * The smallest positive - * number that will be - * entered into the global - * matrix. All smaller - * absolute values will be - * treated as zero and will - * not be assembled. - */ - const double threshold; - - }; - - -/** - * Assemble local matrices into level matrices without using - * block structure. - * - * @todo The matrix structures needed for assembling level matrices - * with local refinement and continuous elements are missing. - * - * @ingroup MeshWorker - * @author Guido Kanschat, 2009 - */ - template - class MGMatrixSimple - { - public: - /** - * Constructor, initializing - * the #threshold, which - * limits how small numbers - * may be to be entered into - * the matrix. - */ - MGMatrixSimple(double threshold = 1.e-12); - - /** - * Store the result matrix - * for later assembling. - */ - void initialize(MGLevelObject& m); - - /** - * Initialize the multilevel - * constraints. - */ - void initialize(const MGConstrainedDoFs& mg_constrained_dofs); - - /** - * Initialize the matrices - * #flux_up and #flux_down - * used for local refinement - * with discontinuous - * Galerkin methods. - */ - void initialize_fluxes(MGLevelObject& flux_up, - MGLevelObject& flux_down); - - /** - * Initialize the matrices - * #interface_in and #interface_out - * used for local refinement - * with continuous - * Galerkin methods. - */ - - void initialize_interfaces(MGLevelObject& interface_in, - MGLevelObject& interface_out); - /** - * Initialize the local data - * in the - * DoFInfo - * object used later for - * assembling. - * - * The info object refers to - * a cell if - * !face, or - * else to an interior or - * boundary face. - */ - template - void initialize_info(DOFINFO& info, bool face) const; - - /** - * Assemble the matrix - * DoFInfo::M1[0] - * into the global matrix. - */ - template - void assemble(const DOFINFO& info); - - /** - * Assemble both local - * matrices in the info - * objects into the global - * matrices. - */ - template - void assemble(const DOFINFO& info1, - const DOFINFO& info2); - private: - /** - * Assemble a single matrix - * into a global matrix. - */ - void assemble(MATRIX& G, - const FullMatrix& M, - const std::vector& i1, - const std::vector& i2); - - /** - * Assemble a single matrix - * into a global matrix. - */ - void assemble(MATRIX& G, - const FullMatrix& M, - const std::vector& i1, - const std::vector& i2, - const unsigned int level); - - /** - * Assemble a single matrix - * into a global matrix. - */ - - void assemble_up(MATRIX& G, - const FullMatrix& M, - const std::vector& i1, - const std::vector& i2, - const unsigned int level = -1); - /** - * Assemble a single matrix - * into a global matrix. - */ - - void assemble_down(MATRIX& G, - const FullMatrix& M, - const std::vector& i1, - const std::vector& i2, - const unsigned int level = -1); - - /** - * Assemble a single matrix - * into a global matrix. - */ - - void assemble_in(MATRIX& G, - const FullMatrix& M, - const std::vector& i1, - const std::vector& i2, - const unsigned int level = -1); - - /** - * Assemble a single matrix - * into a global matrix. - */ - - void assemble_out(MATRIX& G, - const FullMatrix& M, - const std::vector& i1, - const std::vector& i2, - const unsigned int level = -1); - - /** - * The global matrix being - * assembled. + * + * The info object refers to + * a cell if + * !face, or + * else to an interior or + * boundary face. */ - SmartPointer,MGMatrixSimple > matrix; + template + void initialize_info(DOFINFO& info, bool face) const; + /** - * The matrix used for face - * flux terms across the - * refinement edge, coupling - * coarse to fine. + * Assemble the local residuals + * into the global residuals. */ - SmartPointer,MGMatrixSimple > flux_up; + template + void assemble(const DOFINFO& info); /** - * The matrix used for face - * flux terms across the - * refinement edge, coupling - * fine to coarse. + * Assemble both local residuals + * into the global residuals. */ - SmartPointer,MGMatrixSimple > flux_down; - + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); + private: /** - * The matrix used for face - * contributions for continuous - * elements across the - * refinement edge, coupling - * coarse to fine. + * Assemble a single local + * residual into the global. */ - SmartPointer,MGMatrixSimple > interface_in; + void assemble(VECTOR& global, + const BlockVector& local, + const std::vector& dof); /** - * The matrix used for face - * contributions for continuous - * elements across the - * refinement edge, coupling - * fine to coarse. + * The global matrices, + * stored as a vector of + * pointers. */ - SmartPointer,MGMatrixSimple > interface_out; + NamedData > > residuals; + + /** + * A pointer to the object containing the block structure. + */ + SmartPointer > block_info; /** * A pointer to the object containing constraints. */ - SmartPointer > mg_constrained_dofs; - - /** - * The smallest positive - * number that will be - * entered into the global - * matrix. All smaller - * absolute values will be - * treated as zero and will - * not be assembled. - */ - const double threshold; - - }; - - + SmartPointer > constraints; + }; + + /** * A helper class assembling local matrices into global matrices. * @@ -947,318 +414,62 @@ namespace MeshWorker const std::vector& dof1, const std::vector& dof2, unsigned int level1, - unsigned int level2, - bool transpose = false); - - /** - * The level matrices, - * stored as a vector of - * pointers. - */ - MatrixPtrVectorPtr matrices; - - /** - * The flux matrix between - * the fine and the coarse - * level at refinement edges. - */ - MatrixPtrVectorPtr flux_down; - - /** - * The flux matrix between - * the coarse and the fine - * level at refinement edges. - */ - MatrixPtrVectorPtr flux_up; - - /** - * The interface matrix between - * the fine and the coarse - * level at refinement edges. - */ - MatrixPtrVectorPtr interface_out; - - /** - * The interface matrix between - * the coarse and the fine - * level at refinement edges. - */ - MatrixPtrVectorPtr interface_in; - - /** - * A pointer to the object containing the block structure. - */ - SmartPointer > block_info; - - /** - * The smallest positive - * number that will be - * entered into the global - * matrix. All smaller - * absolute values will be - * treated as zero and will - * not be assembled. - */ - const double threshold; - - }; - -/** - * Assemble a simple matrix and a simple right hand side at once. We - * use a combination of MatrixSimple and ResidualSimple to achieve - * this. Cell and face operators should fill the matrix and vector - * objects in LocalResults and this class will assemble - * them into matrix and vector objects. - * - * @ingroup MeshWorker - * @author Guido Kanschat, 2009 - */ - template - class SystemSimple : - private MatrixSimple, - private ResidualSimple - { - public: - /** - * Constructor setting the - * threshold value in - * MatrixSimple. - */ - SystemSimple(double threshold = 1.e-12); - - /** - * Store the two objects data - * is assembled into. - */ - void initialize(MATRIX& m, VECTOR& rhs); - - /** - * Initialize the local data - * in the - * DoFInfo - * object used later for - * assembling. - * - * The info object refers to - * a cell if - * !face, or - * else to an interior or - * boundary face. - */ - template - void initialize_info(DOFINFO& info, bool face) const; - - /** - * Assemble the matrix - * DoFInfo::M1[0] - * into the global matrix. - */ - template - void assemble(const DOFINFO& info); - - /** - * Assemble both local - * matrices in the info - * objects into the global - * matrix. - */ - template - void assemble(const DOFINFO& info1, - const DOFINFO& info2); - }; - - -//----------------------------------------------------------------------// - - template - inline void - Functional::initialize(unsigned int n) - { - results.resize(n); - std::fill(results.begin(), results.end(), 0.); - } - - - template - template - inline void - Functional::initialize_info(DOFINFO& info, bool) - { - info.initialize_numbers(results.size()); - } - - - template - template - inline void - Functional::assemble(const DOFINFO& info) - { - for (unsigned int i=0;i - template - inline void - Functional::assemble(const DOFINFO& info1, - const DOFINFO& info2) - { - for (unsigned int i=0;i - inline number - Functional::operator() (unsigned int i) const - { - AssertIndexRange(i, results.size()); - return results[i]; - } - -//----------------------------------------------------------------------// - - template - inline void - CellsAndFaces::initialize(DataVectors& r, bool sep) - { - Assert(r.name(0) == "cells", typename DataVectors::ExcNameMismatch(0, "cells")); - if (sep) - { - Assert(r.name(1) == "faces", typename DataVectors::ExcNameMismatch(1, "faces")); - AssertDimension(r(0)->n_blocks(), r(1)->n_blocks()); - } - - results = r; - separate_faces = sep; - } - - - template - template - inline void - CellsAndFaces::initialize_info(DOFINFO& info, bool) const - { - info.initialize_numbers(results(0)->n_blocks()); - } - - - template - template - inline void - CellsAndFaces::assemble(const DOFINFO& info) - { - for (unsigned int i=0;iblock(i)(info.face->user_index()) += info.value(i); - else - results(0)->block(i)(info.cell->user_index()) += info.value(i); - } - } - - - template - template - inline void - CellsAndFaces::assemble(const DOFINFO& info1, - const DOFINFO& info2) - { - for (unsigned int i=0;iblock(i)(info1.face->user_index()) += J; - if (info2.face != info1.face) - results(1)->block(i)(info2.face->user_index()) += J; - } - else - { - results(0)->block(i)(info1.cell->user_index()) += .5*info1.value(i); - results(0)->block(i)(info2.cell->user_index()) += .5*info2.value(i); - } - } - } - - - -//----------------------------------------------------------------------// - - template - inline void - ResidualSimple::initialize(NamedData& results) - { - residuals = results; - } - - template - inline void - ResidualSimple::initialize(const ConstraintMatrix& c) - { - constraints = &c; - } - - - template - template - inline void - ResidualSimple::initialize_info(DOFINFO& info, bool) const - { - info.initialize_vectors(residuals.size()); - } - - - template - template - inline void - ResidualSimple::assemble(const DOFINFO& info) - { - for (unsigned int k=0;kdistribute_local_to_global( - info.vector(k).block(0), info.indices, (*residuals(k))); - } - } + unsigned int level2, + bool transpose = false); + + /** + * The level matrices, + * stored as a vector of + * pointers. + */ + MatrixPtrVectorPtr matrices; + + /** + * The flux matrix between + * the fine and the coarse + * level at refinement edges. + */ + MatrixPtrVectorPtr flux_down; + + /** + * The flux matrix between + * the coarse and the fine + * level at refinement edges. + */ + MatrixPtrVectorPtr flux_up; - - template - template - inline void - ResidualSimple::assemble(const DOFINFO& info1, - const DOFINFO& info2) - { - for (unsigned int k=0;kdistribute_local_to_global( - info1.vector(k).block(0), info1.indices, (*residuals(k))); - constraints->distribute_local_to_global( - info2.vector(k).block(0), info2.indices, (*residuals(k))); - } - } - } + /** + * The interface matrix between + * the fine and the coarse + * level at refinement edges. + */ + MatrixPtrVectorPtr interface_out; + + /** + * The interface matrix between + * the coarse and the fine + * level at refinement edges. + */ + MatrixPtrVectorPtr interface_in; + + /** + * A pointer to the object containing the block structure. + */ + SmartPointer > block_info; + + /** + * The smallest positive + * number that will be + * entered into the global + * matrix. All smaller + * absolute values will be + * treated as zero and will + * not be assembled. + */ + const double threshold; + + }; - //----------------------------------------------------------------------// template @@ -1343,421 +554,6 @@ namespace MeshWorker } -//----------------------------------------------------------------------// - - template - inline - MatrixSimple::MatrixSimple(double threshold) - : - threshold(threshold) - {} - - - template - inline void - MatrixSimple::initialize(MATRIX& m) - { - matrix = &m; - } - - - template - inline void - MatrixSimple::initialize(const ConstraintMatrix& c) - { - constraints = &c; - } - - - template - template - inline void - MatrixSimple::initialize_info(DOFINFO& info, bool face) const - { - info.initialize_matrices(1, face); - } - - - - template - inline void - MatrixSimple::assemble(const FullMatrix& M, - const std::vector& i1, - const std::vector& i2) - { - AssertDimension(M.m(), i1.size()); - AssertDimension(M.n(), i2.size()); - - if(constraints == 0) - { - for (unsigned int j=0; j= threshold) - matrix->add(i1[j], i2[k], M(j,k)); - } - else - constraints->distribute_local_to_global( - M, i1, i2, *matrix); - } - - - template - template - inline void - MatrixSimple::assemble(const DOFINFO& info) - { - assemble(info.matrix(0,false).matrix, info.indices, info.indices); - } - - - template - template - inline void - MatrixSimple::assemble(const DOFINFO& info1, - const DOFINFO& info2) - { - assemble(info1.matrix(0,false).matrix, info1.indices, info1.indices); - assemble(info1.matrix(0,true).matrix, info1.indices, info2.indices); - assemble(info2.matrix(0,false).matrix, info2.indices, info2.indices); - assemble(info2.matrix(0,true).matrix, info2.indices, info1.indices); - } - - -//----------------------------------------------------------------------// - - template - inline - MGMatrixSimple::MGMatrixSimple(double threshold) - : - threshold(threshold) - {} - - - template - inline void - MGMatrixSimple::initialize(MGLevelObject& m) - { - matrix = &m; - } - - template - inline void - MGMatrixSimple::initialize(const MGConstrainedDoFs& c) - { - mg_constrained_dofs = &c; - } - - template - inline void - MGMatrixSimple::initialize_fluxes( - MGLevelObject& up, MGLevelObject& down) - { - flux_up = &up; - flux_down = &down; - } - - - template - inline void - MGMatrixSimple::initialize_interfaces( - MGLevelObject& in, MGLevelObject& out) - { - interface_in = ∈ - interface_out = &out; - } - - - template - template - inline void - MGMatrixSimple::initialize_info(DOFINFO& info, bool face) const - { - info.initialize_matrices(1, face); - } - - - template - inline void - MGMatrixSimple::assemble( - MATRIX& G, - const FullMatrix& M, - const std::vector& i1, - const std::vector& i2) - { - AssertDimension(M.m(), i1.size()); - AssertDimension(M.n(), i2.size()); - - if(mg_constrained_dofs == 0) - { - for (unsigned int j=0; j= threshold) - G.add(i1[j], i2[k], M(j,k)); - } - else - { - for (unsigned int j=0; j= threshold) - { - if(!mg_constrained_dofs->continuity_across_refinement_edges()) - G.add(i1[j], i2[k], M(j,k)); - } - } - } - - - template - inline void - MGMatrixSimple::assemble( - MATRIX& G, - const FullMatrix& M, - const std::vector& i1, - const std::vector& i2, - const unsigned int level) - { - AssertDimension(M.m(), i1.size()); - AssertDimension(M.n(), i2.size()); - - if(mg_constrained_dofs == 0) - { - for (unsigned int j=0; j= threshold) - G.add(i1[j], i2[k], M(j,k)); - } - else - { - for (unsigned int j=0; j= threshold) - if (!mg_constrained_dofs->at_refinement_edge(level, i1[j]) && - !mg_constrained_dofs->at_refinement_edge(level, i2[k])) - { - if (mg_constrained_dofs->set_boundary_values()) - { - if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) && - !mg_constrained_dofs->is_boundary_index(level, i2[k])) - || - (mg_constrained_dofs->is_boundary_index(level, i1[j]) && - mg_constrained_dofs->is_boundary_index(level, i2[k]) && - i1[j] == i2[k])) - G.add(i1[j], i2[k], M(j,k)); - } - else - G.add(i1[j], i2[k], M(j,k)); - } - } - } - - - template - inline void - MGMatrixSimple::assemble_up( - MATRIX& G, - const FullMatrix& M, - const std::vector& i1, - const std::vector& i2, - const unsigned int level) - { - AssertDimension(M.n(), i1.size()); - AssertDimension(M.m(), i2.size()); - - if(mg_constrained_dofs == 0) - { - for (unsigned int j=0; j= threshold) - G.add(i1[j], i2[k], M(k,j)); - } - else - { - for (unsigned int j=0; j= threshold) - if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) && - !mg_constrained_dofs->at_refinement_edge(level, i2[k])) - G.add(i1[j], i2[k], M(k,j)); - } - } - - template - inline void - MGMatrixSimple::assemble_down( - MATRIX& G, - const FullMatrix& M, - const std::vector& i1, - const std::vector& i2, - const unsigned int level) - { - AssertDimension(M.m(), i1.size()); - AssertDimension(M.n(), i2.size()); - - if(mg_constrained_dofs == 0) - { - for (unsigned int j=0; j= threshold) - G.add(i1[j], i2[k], M(j,k)); - } - else - { - for (unsigned int j=0; j= threshold) - if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) && - !mg_constrained_dofs->at_refinement_edge(level, i2[k])) - G.add(i1[j], i2[k], M(j,k)); - } - } - - template - inline void - MGMatrixSimple::assemble_in( - MATRIX& G, - const FullMatrix& M, - const std::vector& i1, - const std::vector& i2, - const unsigned int level) - { - AssertDimension(M.m(), i1.size()); - AssertDimension(M.n(), i2.size()); - - if(mg_constrained_dofs == 0) - { - for (unsigned int j=0; j= threshold) - G.add(i1[j], i2[k], M(j,k)); - } - else - { - for (unsigned int j=0; j= threshold) - if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) && - !mg_constrained_dofs->at_refinement_edge(level, i2[k])) - { - if (mg_constrained_dofs->set_boundary_values()) - { - if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && - !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])) - || - (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && - mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) && - i1[j] == i2[k])) - G.add(i1[j], i2[k], M(j,k)); - } - else - G.add(i1[j], i2[k], M(j,k)); - } - } - } - - template - inline void - MGMatrixSimple::assemble_out( - MATRIX& G, - const FullMatrix& M, - const std::vector& i1, - const std::vector& i2, - const unsigned int level) - { - AssertDimension(M.n(), i1.size()); - AssertDimension(M.m(), i2.size()); - - if(mg_constrained_dofs == 0) - { - for (unsigned int j=0; j= threshold) - G.add(i1[j], i2[k], M(k,j)); - } - else - { - for (unsigned int j=0; j= threshold) - if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) && - !mg_constrained_dofs->at_refinement_edge(level, i2[k])) - { - if (mg_constrained_dofs->set_boundary_values()) - { - if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && - !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])) - || - (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && - mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) && - i1[j] == i2[k])) - G.add(i1[j], i2[k], M(k,j)); - } - else - G.add(i1[j], i2[k], M(k,j)); - } - } - } - - - template - template - inline void - MGMatrixSimple::assemble(const DOFINFO& info) - { - const unsigned int level = info.cell->level(); - assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices, level); - - if(mg_constrained_dofs != 0) - { - assemble_in((*interface_in)[level], info.matrix(0,false).matrix, info.indices, info.indices, level); - assemble_out((*interface_out)[level],info.matrix(0,false).matrix, info.indices, info.indices, level); - } - } - - - template - template - inline void - MGMatrixSimple::assemble(const DOFINFO& info1, - const DOFINFO& info2) - { - const unsigned int level1 = info1.cell->level(); - const unsigned int level2 = info2.cell->level(); - - if (level1 == level2) - { - if(mg_constrained_dofs == 0) - { - assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices); - assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices); - assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices); - assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices); - } - else - { - assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1); - assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1); - assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1); - assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1); - } - } - else - { - Assert(level1 > level2, ExcInternalError()); - // Do not add info2.M1, - // which is done by - // the coarser cell - assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices); - //assemble_transpose((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices); - //assemble((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices); - if(level1>0) - { - assemble_up((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices, level1); - assemble_down((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1); - } - } - } - - //----------------------------------------------------------------------// template @@ -2052,61 +848,6 @@ namespace MeshWorker } } -//----------------------------------------------------------------------// - - template - SystemSimple::SystemSimple(double t) - : - MatrixSimple(t) - {} - - - template - inline void - SystemSimple::initialize(MATRIX& m, VECTOR& rhs) - { - NamedData data; - VECTOR* p = &rhs; - data.add(p, "right hand side"); - - MatrixSimple::initialize(m); - ResidualSimple::initialize(data); - } - - - template - template - inline void - SystemSimple::initialize_info(DOFINFO& info, - bool face) const - { - MatrixSimple::initialize_info(info, face); - ResidualSimple::initialize_info(info, face); - } - - - template - template - inline void - SystemSimple::assemble(const DOFINFO& info) - { - MatrixSimple::assemble(info); - ResidualSimple::assemble(info); - } - - - template - template - inline void - SystemSimple::assemble(const DOFINFO& info1, - const DOFINFO& info2) - { - MatrixSimple::assemble(info1, info2); - ResidualSimple::assemble(info1, info2); - } - } -} - DEAL_II_NAMESPACE_CLOSE #endif diff --git a/deal.II/include/deal.II/meshworker/functional.h b/deal.II/include/deal.II/meshworker/functional.h new file mode 100644 index 0000000000..3a7417436c --- /dev/null +++ b/deal.II/include/deal.II/meshworker/functional.h @@ -0,0 +1,311 @@ +//--------------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2010, 2011 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__mesh_worker_functional_h +#define __deal2__mesh_worker_functional_h + +#include +#include +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace MeshWorker +{ + namespace Assembler + { +/** + * The class assembling local contributions to a functional into + * global functionals. + * + * + * + * @ingroup MeshWorker + * @author Guido Kanschat, 2009 + */ + template + class Functional + { + public: + /** + * Initialize local data to + * store functionals. The + * number n is the + * number of functionals to + * be computed. + */ + void initialize(unsigned int n); + /** + * Initialize the local data + * in the DoFInfo object used + * later for assembling. + * + * The info object refers to + * a cell if + * !face, or + * else to an interior or + * boundary face. + */ + template + void initialize_info(DOFINFO& info, bool face); + + /** + * Assemble the local values + * into the global vectors. + */ + template + void assemble(const DOFINFO& info); + + /** + * Assemble both local values + * into the global vectors. + */ + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); + + /** + * The value of the ith entry + * in #results. + */ + number operator() (unsigned int i) const; + private: + /** + * The values into which the + * results are added. + */ + std::vector results; + }; + +/** + * Compute cell and face contributions of one or several functionals, + * typically for error estimates. + * + * @ingroup MeshWorker + * @author Guido Kanschat, 2009 + */ + template + class CellsAndFaces + { + public: + /** + * The data type for + * communicating the cell and + * face vectors. + */ + typedef NamedData*> DataVectors; + + /** + * The initialization + * function, specifying the + * #results vectors and + * whether face data should + * be collected separately. + * + * #results should contain + * two block vectors named + * "cells" and "faces" (the + * latter only if + * #separate_faces is + * true). In each of the two, + * each block should have + * equal size and be large + * enough to accomodate all + * user indices set in the + * cells and faces covered by + * the loop it is used + * in. Typically, for + * estimators, this is + * Triangulation::n_active_cells() + * and + * Triangulation::n_faces(), + * respectively. + * + * The use of BlockVector may + * seem cumbersome, but it + * allows us to assemble + * several functionals at the + * same time, one in each + * block. The typical + * situation for error + * estimate is just having a + * single block in each vector. + */ + void initialize(DataVectors& results, + bool separate_faces = true); + /** + * Initialize the local data + * in the + * DoFInfo + * object used later for + * assembling. + * + * The info object refers to + * a cell if + * !face, or + * else to an interior or + * boundary face. + */ + template + void initialize_info(DOFINFO& info, bool face) const; + + /** + * Assemble the local values + * into the global vectors. + */ + template + void assemble(const DOFINFO& info); + + /** + * Assemble both local values + * into the global vectors. + */ + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); + + /** + * The value of the ith entry + * in #results. + */ + number operator() (unsigned int i) const; + private: + DataVectors results; + bool separate_faces; + }; +//----------------------------------------------------------------------// + + template + inline void + Functional::initialize(unsigned int n) + { + results.resize(n); + std::fill(results.begin(), results.end(), 0.); + } + + + template + template + inline void + Functional::initialize_info(DOFINFO& info, bool) + { + info.initialize_numbers(results.size()); + } + + + template + template + inline void + Functional::assemble(const DOFINFO& info) + { + for (unsigned int i=0;i + template + inline void + Functional::assemble(const DOFINFO& info1, + const DOFINFO& info2) + { + for (unsigned int i=0;i + inline number + Functional::operator() (unsigned int i) const + { + AssertIndexRange(i, results.size()); + return results[i]; + } + +//----------------------------------------------------------------------// + + template + inline void + CellsAndFaces::initialize(DataVectors& r, bool sep) + { + Assert(r.name(0) == "cells", typename DataVectors::ExcNameMismatch(0, "cells")); + if (sep) + { + Assert(r.name(1) == "faces", typename DataVectors::ExcNameMismatch(1, "faces")); + AssertDimension(r(0)->n_blocks(), r(1)->n_blocks()); + } + + results = r; + separate_faces = sep; + } + + + template + template + inline void + CellsAndFaces::initialize_info(DOFINFO& info, bool) const + { + info.initialize_numbers(results(0)->n_blocks()); + } + + + template + template + inline void + CellsAndFaces::assemble(const DOFINFO& info) + { + for (unsigned int i=0;iblock(i)(info.face->user_index()) += info.value(i); + else + results(0)->block(i)(info.cell->user_index()) += info.value(i); + } + } + + + template + template + inline void + CellsAndFaces::assemble(const DOFINFO& info1, + const DOFINFO& info2) + { + for (unsigned int i=0;iblock(i)(info1.face->user_index()) += J; + if (info2.face != info1.face) + results(1)->block(i)(info2.face->user_index()) += J; + } + else + { + results(0)->block(i)(info1.cell->user_index()) += .5*info1.value(i); + results(0)->block(i)(info2.cell->user_index()) += .5*info2.value(i); + } + } + } + } +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/include/deal.II/meshworker/simple.h b/deal.II/include/deal.II/meshworker/simple.h new file mode 100644 index 0000000000..81970cce62 --- /dev/null +++ b/deal.II/include/deal.II/meshworker/simple.h @@ -0,0 +1,1021 @@ +//--------------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2010, 2011 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__mesh_worker_assembler_h +#define __deal2__mesh_worker_assembler_h + +#include +#include +#include +#include +#include +#include +#include + +/** + * @file + * @brief MeshWorker::Assember objects for systems without block structure + * + * The header containing the classes + * MeshWorker::Assember::MatrixSimple, + * MeshWorker::Assember::MGMatrixSimple, + * MeshWorker::Assember::ResidualSimple, and + * MeshWorker::Assember::SystemSimple. + */ + +DEAL_II_NAMESPACE_OPEN + +namespace MeshWorker +{ + namespace Assembler + { +/** + * Assemble residuals without block structure. + * + * The data structure for this Assembler class is a simple vector on + * each cell with entries from zero to + * FiniteElementData::dofs_per_cell and a simple global vector with + * entries numbered from zero to DoFHandler::n_dofs(). No BlockInfo is + * required and the global vector may be any type of vector having + * element access through operator() (unsigned int) + * + * @ingroup MeshWorker + * @author Guido Kanschat, 2009 + */ + template + class ResidualSimple + { + public: + void initialize(NamedData& results); + /** + * Initialize the constraints. + */ + void initialize(const ConstraintMatrix& constraints); + /** + * Initialize the local data + * in the DoFInfo object used + * later for assembling. + * + * The info object refers to + * a cell if + * !face, or + * else to an interior or + * boundary face. + */ + template + void initialize_info(DOFINFO& info, bool face) const; + + /** + * Assemble the local residuals + * into the global residuals. + */ + template + void assemble(const DOFINFO& info); + + /** + * Assemble both local residuals + * into the global residuals. + */ + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); + private: + /** + * The global residal vectors + * filled by assemble(). + */ + NamedData > > residuals; + /** + * A pointer to the object containing constraints. + */ + SmartPointer > constraints; + }; + +/** + * Assemble local matrices into a single global matrix without using + * block structure. + * + * After being initialized with a SparseMatrix object (or another + * matrix offering the same functionality as SparseMatrix::add()), + * this class can be used in a MeshWorker::loop() to assemble the cell + * and face matrices into the global matrix. + * + * @todo On locally refined meshes, a ConstraintMatrix should be used + * to automatically eliminate hanging nodes. + * + * @ingroup MeshWorker + * @author Guido Kanschat, 2009 + */ + template + class MatrixSimple + { + public: + /** + * Constructor, initializing + * the #threshold, which + * limits how small numbers + * may be to be entered into + * the matrix. + */ + MatrixSimple(double threshold = 1.e-12); + + /** + * Store the result matrix + * for later assembling. + */ + void initialize(MATRIX& m); + /** + * Initialize the constraints. + */ + void initialize(const ConstraintMatrix& constraints); + /** + * Initialize the local data + * in the + * DoFInfo + * object used later for + * assembling. + * + * The info object refers to + * a cell if + * !face, or + * else to an interior or + * boundary face. + */ + template + void initialize_info(DOFINFO& info, bool face) const; + + /** + * Assemble the matrix + * DoFInfo::M1[0] + * into the global matrix. + */ + template + void assemble(const DOFINFO& info); + + /** + * Assemble both local + * matrices in the info + * objects into the global + * matrix. + */ + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); + private: + /** + * Assemble a single matrix + * into #matrix. + */ + void assemble(const FullMatrix& M, + const std::vector& i1, + const std::vector& i2); + + /** + * The global matrix being + * assembled. + */ + SmartPointer > matrix; + /** + * A pointer to the object containing constraints. + */ + SmartPointer > constraints; + + /** + * The smallest positive + * number that will be + * entered into the global + * matrix. All smaller + * absolute values will be + * treated as zero and will + * not be assembled. + */ + const double threshold; + + }; + + +/** + * Assemble local matrices into level matrices without using + * block structure. + * + * @todo The matrix structures needed for assembling level matrices + * with local refinement and continuous elements are missing. + * + * @ingroup MeshWorker + * @author Guido Kanschat, 2009 + */ + template + class MGMatrixSimple + { + public: + /** + * Constructor, initializing + * the #threshold, which + * limits how small numbers + * may be to be entered into + * the matrix. + */ + MGMatrixSimple(double threshold = 1.e-12); + + /** + * Store the result matrix + * for later assembling. + */ + void initialize(MGLevelObject& m); + + /** + * Initialize the multilevel + * constraints. + */ + void initialize(const MGConstrainedDoFs& mg_constrained_dofs); + + /** + * Initialize the matrices + * #flux_up and #flux_down + * used for local refinement + * with discontinuous + * Galerkin methods. + */ + void initialize_fluxes(MGLevelObject& flux_up, + MGLevelObject& flux_down); + + /** + * Initialize the matrices + * #interface_in and #interface_out + * used for local refinement + * with continuous + * Galerkin methods. + */ + + void initialize_interfaces(MGLevelObject& interface_in, + MGLevelObject& interface_out); + /** + * Initialize the local data + * in the + * DoFInfo + * object used later for + * assembling. + * + * The info object refers to + * a cell if + * !face, or + * else to an interior or + * boundary face. + */ + template + void initialize_info(DOFINFO& info, bool face) const; + + /** + * Assemble the matrix + * DoFInfo::M1[0] + * into the global matrix. + */ + template + void assemble(const DOFINFO& info); + + /** + * Assemble both local + * matrices in the info + * objects into the global + * matrices. + */ + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); + private: + /** + * Assemble a single matrix + * into a global matrix. + */ + void assemble(MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2); + + /** + * Assemble a single matrix + * into a global matrix. + */ + void assemble(MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2, + const unsigned int level); + + /** + * Assemble a single matrix + * into a global matrix. + */ + + void assemble_up(MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2, + const unsigned int level = -1); + /** + * Assemble a single matrix + * into a global matrix. + */ + + void assemble_down(MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2, + const unsigned int level = -1); + + /** + * Assemble a single matrix + * into a global matrix. + */ + + void assemble_in(MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2, + const unsigned int level = -1); + + /** + * Assemble a single matrix + * into a global matrix. + */ + + void assemble_out(MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2, + const unsigned int level = -1); + + /** + * The global matrix being + * assembled. + */ + SmartPointer,MGMatrixSimple > matrix; + + /** + * The matrix used for face + * flux terms across the + * refinement edge, coupling + * coarse to fine. + */ + SmartPointer,MGMatrixSimple > flux_up; + + /** + * The matrix used for face + * flux terms across the + * refinement edge, coupling + * fine to coarse. + */ + SmartPointer,MGMatrixSimple > flux_down; + + /** + * The matrix used for face + * contributions for continuous + * elements across the + * refinement edge, coupling + * coarse to fine. + */ + SmartPointer,MGMatrixSimple > interface_in; + + /** + * The matrix used for face + * contributions for continuous + * elements across the + * refinement edge, coupling + * fine to coarse. + */ + SmartPointer,MGMatrixSimple > interface_out; + /** + * A pointer to the object containing constraints. + */ + SmartPointer > mg_constrained_dofs; + + /** + * The smallest positive + * number that will be + * entered into the global + * matrix. All smaller + * absolute values will be + * treated as zero and will + * not be assembled. + */ + const double threshold; + + }; + + +/** + * Assemble a simple matrix and a simple right hand side at once. We + * use a combination of MatrixSimple and ResidualSimple to achieve + * this. Cell and face operators should fill the matrix and vector + * objects in LocalResults and this class will assemble + * them into matrix and vector objects. + * + * @ingroup MeshWorker + * @author Guido Kanschat, 2009 + */ + template + class SystemSimple : + private MatrixSimple, + private ResidualSimple + { + public: + /** + * Constructor setting the + * threshold value in + * MatrixSimple. + */ + SystemSimple(double threshold = 1.e-12); + + /** + * Store the two objects data + * is assembled into. + */ + void initialize(MATRIX& m, VECTOR& rhs); + + /** + * Initialize the local data + * in the + * DoFInfo + * object used later for + * assembling. + * + * The info object refers to + * a cell if + * !face, or + * else to an interior or + * boundary face. + */ + template + void initialize_info(DOFINFO& info, bool face) const; + + /** + * Assemble the matrix + * DoFInfo::M1[0] + * into the global matrix. + */ + template + void assemble(const DOFINFO& info); + + /** + * Assemble both local + * matrices in the info + * objects into the global + * matrix. + */ + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); + }; + + +//----------------------------------------------------------------------// + + template + inline void + ResidualSimple::initialize(NamedData& results) + { + residuals = results; + } + + template + inline void + ResidualSimple::initialize(const ConstraintMatrix& c) + { + constraints = &c; + } + + + template + template + inline void + ResidualSimple::initialize_info(DOFINFO& info, bool) const + { + info.initialize_vectors(residuals.size()); + } + + + template + template + inline void + ResidualSimple::assemble(const DOFINFO& info) + { + for (unsigned int k=0;kdistribute_local_to_global( + info.vector(k).block(0), info.indices, (*residuals(k))); + } + } + + + template + template + inline void + ResidualSimple::assemble(const DOFINFO& info1, + const DOFINFO& info2) + { + for (unsigned int k=0;kdistribute_local_to_global( + info1.vector(k).block(0), info1.indices, (*residuals(k))); + constraints->distribute_local_to_global( + info2.vector(k).block(0), info2.indices, (*residuals(k))); + } + } + } + + +//----------------------------------------------------------------------// + + template + inline + MatrixSimple::MatrixSimple(double threshold) + : + threshold(threshold) + {} + + + template + inline void + MatrixSimple::initialize(MATRIX& m) + { + matrix = &m; + } + + + template + inline void + MatrixSimple::initialize(const ConstraintMatrix& c) + { + constraints = &c; + } + + + template + template + inline void + MatrixSimple::initialize_info(DOFINFO& info, bool face) const + { + info.initialize_matrices(1, face); + } + + + + template + inline void + MatrixSimple::assemble(const FullMatrix& M, + const std::vector& i1, + const std::vector& i2) + { + AssertDimension(M.m(), i1.size()); + AssertDimension(M.n(), i2.size()); + + if(constraints == 0) + { + for (unsigned int j=0; j= threshold) + matrix->add(i1[j], i2[k], M(j,k)); + } + else + constraints->distribute_local_to_global( + M, i1, i2, *matrix); + } + + + template + template + inline void + MatrixSimple::assemble(const DOFINFO& info) + { + assemble(info.matrix(0,false).matrix, info.indices, info.indices); + } + + + template + template + inline void + MatrixSimple::assemble(const DOFINFO& info1, + const DOFINFO& info2) + { + assemble(info1.matrix(0,false).matrix, info1.indices, info1.indices); + assemble(info1.matrix(0,true).matrix, info1.indices, info2.indices); + assemble(info2.matrix(0,false).matrix, info2.indices, info2.indices); + assemble(info2.matrix(0,true).matrix, info2.indices, info1.indices); + } + + +//----------------------------------------------------------------------// + + template + inline + MGMatrixSimple::MGMatrixSimple(double threshold) + : + threshold(threshold) + {} + + + template + inline void + MGMatrixSimple::initialize(MGLevelObject& m) + { + matrix = &m; + } + + template + inline void + MGMatrixSimple::initialize(const MGConstrainedDoFs& c) + { + mg_constrained_dofs = &c; + } + + template + inline void + MGMatrixSimple::initialize_fluxes( + MGLevelObject& up, MGLevelObject& down) + { + flux_up = &up; + flux_down = &down; + } + + + template + inline void + MGMatrixSimple::initialize_interfaces( + MGLevelObject& in, MGLevelObject& out) + { + interface_in = ∈ + interface_out = &out; + } + + + template + template + inline void + MGMatrixSimple::initialize_info(DOFINFO& info, bool face) const + { + info.initialize_matrices(1, face); + } + + + template + inline void + MGMatrixSimple::assemble( + MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2) + { + AssertDimension(M.m(), i1.size()); + AssertDimension(M.n(), i2.size()); + + if(mg_constrained_dofs == 0) + { + for (unsigned int j=0; j= threshold) + G.add(i1[j], i2[k], M(j,k)); + } + else + { + for (unsigned int j=0; j= threshold) + { + if(!mg_constrained_dofs->continuity_across_refinement_edges()) + G.add(i1[j], i2[k], M(j,k)); + } + } + } + + + template + inline void + MGMatrixSimple::assemble( + MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2, + const unsigned int level) + { + AssertDimension(M.m(), i1.size()); + AssertDimension(M.n(), i2.size()); + + if(mg_constrained_dofs == 0) + { + for (unsigned int j=0; j= threshold) + G.add(i1[j], i2[k], M(j,k)); + } + else + { + for (unsigned int j=0; j= threshold) + if (!mg_constrained_dofs->at_refinement_edge(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge(level, i2[k])) + { + if (mg_constrained_dofs->set_boundary_values()) + { + if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) && + !mg_constrained_dofs->is_boundary_index(level, i2[k])) + || + (mg_constrained_dofs->is_boundary_index(level, i1[j]) && + mg_constrained_dofs->is_boundary_index(level, i2[k]) && + i1[j] == i2[k])) + G.add(i1[j], i2[k], M(j,k)); + } + else + G.add(i1[j], i2[k], M(j,k)); + } + } + } + + + template + inline void + MGMatrixSimple::assemble_up( + MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2, + const unsigned int level) + { + AssertDimension(M.n(), i1.size()); + AssertDimension(M.m(), i2.size()); + + if(mg_constrained_dofs == 0) + { + for (unsigned int j=0; j= threshold) + G.add(i1[j], i2[k], M(k,j)); + } + else + { + for (unsigned int j=0; j= threshold) + if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge(level, i2[k])) + G.add(i1[j], i2[k], M(k,j)); + } + } + + template + inline void + MGMatrixSimple::assemble_down( + MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2, + const unsigned int level) + { + AssertDimension(M.m(), i1.size()); + AssertDimension(M.n(), i2.size()); + + if(mg_constrained_dofs == 0) + { + for (unsigned int j=0; j= threshold) + G.add(i1[j], i2[k], M(j,k)); + } + else + { + for (unsigned int j=0; j= threshold) + if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge(level, i2[k])) + G.add(i1[j], i2[k], M(j,k)); + } + } + + template + inline void + MGMatrixSimple::assemble_in( + MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2, + const unsigned int level) + { + AssertDimension(M.m(), i1.size()); + AssertDimension(M.n(), i2.size()); + + if(mg_constrained_dofs == 0) + { + for (unsigned int j=0; j= threshold) + G.add(i1[j], i2[k], M(j,k)); + } + else + { + for (unsigned int j=0; j= threshold) + if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge(level, i2[k])) + { + if (mg_constrained_dofs->set_boundary_values()) + { + if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])) + || + (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && + mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) && + i1[j] == i2[k])) + G.add(i1[j], i2[k], M(j,k)); + } + else + G.add(i1[j], i2[k], M(j,k)); + } + } + } + + template + inline void + MGMatrixSimple::assemble_out( + MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2, + const unsigned int level) + { + AssertDimension(M.n(), i1.size()); + AssertDimension(M.m(), i2.size()); + + if(mg_constrained_dofs == 0) + { + for (unsigned int j=0; j= threshold) + G.add(i1[j], i2[k], M(k,j)); + } + else + { + for (unsigned int j=0; j= threshold) + if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge(level, i2[k])) + { + if (mg_constrained_dofs->set_boundary_values()) + { + if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && + !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])) + || + (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) && + mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) && + i1[j] == i2[k])) + G.add(i1[j], i2[k], M(k,j)); + } + else + G.add(i1[j], i2[k], M(k,j)); + } + } + } + + + template + template + inline void + MGMatrixSimple::assemble(const DOFINFO& info) + { + const unsigned int level = info.cell->level(); + assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices, level); + + if(mg_constrained_dofs != 0) + { + assemble_in((*interface_in)[level], info.matrix(0,false).matrix, info.indices, info.indices, level); + assemble_out((*interface_out)[level],info.matrix(0,false).matrix, info.indices, info.indices, level); + } + } + + + template + template + inline void + MGMatrixSimple::assemble(const DOFINFO& info1, + const DOFINFO& info2) + { + const unsigned int level1 = info1.cell->level(); + const unsigned int level2 = info2.cell->level(); + + if (level1 == level2) + { + if(mg_constrained_dofs == 0) + { + assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices); + assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices); + assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices); + assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices); + } + else + { + assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1); + assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1); + assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1); + assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1); + } + } + else + { + Assert(level1 > level2, ExcInternalError()); + // Do not add info2.M1, + // which is done by + // the coarser cell + assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices); + //assemble_transpose((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices); + //assemble((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices); + if(level1>0) + { + assemble_up((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices, level1); + assemble_down((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1); + } + } + } +//----------------------------------------------------------------------// + + template + SystemSimple::SystemSimple(double t) + : + MatrixSimple(t) + {} + + + template + inline void + SystemSimple::initialize(MATRIX& m, VECTOR& rhs) + { + NamedData data; + VECTOR* p = &rhs; + data.add(p, "right hand side"); + + MatrixSimple::initialize(m); + ResidualSimple::initialize(data); + } + + + template + template + inline void + SystemSimple::initialize_info(DOFINFO& info, + bool face) const + { + MatrixSimple::initialize_info(info, face); + ResidualSimple::initialize_info(info, face); + } + + + template + template + inline void + SystemSimple::assemble(const DOFINFO& info) + { + MatrixSimple::assemble(info); + ResidualSimple::assemble(info); + } + + + template + template + inline void + SystemSimple::assemble(const DOFINFO& info1, + const DOFINFO& info2) + { + MatrixSimple::assemble(info1, info2); + ResidualSimple::assemble(info1, info2); + } + } +} + +DEAL_II_NAMESPACE_CLOSE + +#endif