From: Wolfgang Bangerth Date: Wed, 11 Mar 1998 16:14:20 +0000 (+0000) Subject: Mostly bugfixes in the assembler and fe part. X-Git-Tag: v8.0.0~23195 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=81ec1641ee7e4a369be912f427f6bc1482593de1;p=dealii.git Mostly bugfixes in the assembler and fe part. git-svn-id: https://svn.dealii.org/trunk@56 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/quadrature.h b/deal.II/base/include/base/quadrature.h index 7f06f7c545..3f432333c8 100644 --- a/deal.II/base/include/base/quadrature.h +++ b/deal.II/base/include/base/quadrature.h @@ -41,6 +41,14 @@ class Quadrature { */ double weight (const unsigned int i) const; + /** + * Exception + */ + DeclException2 (ExcInvalidIndex, + int, int, + << "The index " << arg1 + << " is out of range, it should be less than " << arg2); + protected: /** * List of quadrature points. To be filled diff --git a/deal.II/base/source/quadrature.cc b/deal.II/base/source/quadrature.cc new file mode 100644 index 0000000000..187cd73266 --- /dev/null +++ b/deal.II/base/source/quadrature.cc @@ -0,0 +1,34 @@ +/* $Id$ */ + +#include + + +template +Quadrature::Quadrature (const unsigned int n_q) : + n_quadrature_points(n_q), + quadrature_points (n_q, Point()), + weights (n_q, 0) {}; + + + +template +const Point & Quadrature::quad_point (const unsigned int i) const { + Assert (i +double Quadrature::weight (const unsigned int i) const { + Assert (i; +template class Quadrature<2>; diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index 0f21f502fb..44c574b6b6 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -12,8 +12,8 @@ QGauss2<1>::QGauss2 () : for (unsigned int i=0; i (xpts[i])); - weights.push_back (wts[i]); + quadrature_points[i] = Point<1>(xpts[i]); + weights[i] = wts[i]; }; }; @@ -36,8 +36,8 @@ QGauss2x4<1>::QGauss2x4 () : for (unsigned int i=0; i (xpts[i])); - weights.push_back (wts[i]); + quadrature_points[i] = Point<1>(xpts[i]); + weights[i] = wts[i]; }; }; @@ -58,8 +58,8 @@ QGauss4<1>::QGauss4 () : for (unsigned int i=0; i (xpts[i])); - weights.push_back (wts[i]); + quadrature_points[i] = Point<1>(xpts[i]); + weights[i] = wts[i]; }; }; @@ -86,8 +86,8 @@ QGauss8<1>::QGauss8 () : for (unsigned int i=0; i (xpts[i])); - weights.push_back (wts[i]); + quadrature_points[i] = Point<1>(xpts[i]); + weights[i] = wts[i]; }; }; @@ -96,8 +96,8 @@ QGauss8<1>::QGauss8 () : QMidpoint<1>::QMidpoint () : Quadrature<1>(1) { - quadrature_points.push_back (0.5); - weights.push_back (1.0); + quadrature_points[0] = 0.5; + weights[0] = 1.0; }; @@ -110,8 +110,8 @@ QSimpson<1>::QSimpson () : for (unsigned int i=0; i (xpts[i])); - weights.push_back (wts[i]); + quadrature_points[i] = Point<1>(xpts[i]); + weights[i] = wts[i]; }; }; @@ -125,7 +125,26 @@ QTrapez<1>::QTrapez () : for (unsigned int i=0; i (xpts[i])); - weights.push_back (wts[i]); + quadrature_points[i] = Point<1>(xpts[i]); + weights[i] = wts[i]; }; }; + + + + +QGauss4<2>::QGauss4 () : + Quadrature<2> (4) { + static const double xpts[] = { 0.211324865, 0.788675135, + 0.211324865, 0.788675135 }; + static const double ypts[] = { 0.211324865, 0.211324865, + 0.788675135, 0.788675135 }; + static const double wts[] = { 1./4., 1./4., 1./4., 1./4. }; + + for (unsigned int i=0; i(xpts[i], ypts[i]); + weights[i] = wts[i]; + }; + +}; diff --git a/deal.II/deal.II/Attic/examples/Makefile b/deal.II/deal.II/Attic/examples/Makefile index 719e01c5d0..e874cc83e6 100644 --- a/deal.II/deal.II/Attic/examples/Makefile +++ b/deal.II/deal.II/Attic/examples/Makefile @@ -5,8 +5,9 @@ INCPATH = -I../include -I../../lac/include -I../../base/include CXXFLAGS.g= -DDEBUG -g -Wall $(INCPATH) CXXFLAGS = -O3 $(INCPATH) -LIBS.g = -L../lib -lbasic.g -lgrid.g -lfe.g -lnumerics.g -llac.g -LIBS = -L../lib -lbasic -lgrid -lfe -lnumerics +LIBPATH = -L../lib -L../../lac/lib +LIBS.g = $(LIBPATH) -lbasic.g -lgrid.g -lfe.g -lnumerics.g -llac.g +LIBS = $(LIBPATH) -lbasic -lgrid -lfe -lnumerics ifeq ($(shell uname),SunOS) CXXFLAGS := -V2.7.2.3 $(CXXFLAGS) @@ -14,17 +15,22 @@ CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g) endif -cc-files = grid/grid_test.cc +cc-files = grid/grid_test.cc dof/dof_test.cc poisson/poisson.cc o-files = $(cc-files:.cc=.o) h-files = $(wildcard ../include/*.h) -examples: grid/grid_test dof/dof_test +examples: $(cc-files:.cc=) +%.o : %.cc + $(CXX) $(CXXFLAGS) -o $@ -c $< -grid/grid_test: grid/grid_test.o ../lib/libgrid.a +grid/grid_test: grid/grid_test.o $(CXX) $(CXXFLAGS) -o $@ $< $(LIBS) -dof/dof_test: dof/dof_test.o ../lib/libgrid.g.a ../lib/libdeal.g.a +dof/dof_test: dof/dof_test.o + $(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS.g) -lg++ + +poisson/poisson: poisson/poisson.o $(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS.g) -lg++ @@ -52,6 +58,6 @@ clean: #to remake this file upon inclusion at the bottom #of this file. Makefile.dep: $(cc-files) $(h-files) - $(CXX) -I../include -M $(cc-files) > Makefile.dep + $(CXX) $(INCPATH) -M $(cc-files) > Makefile.dep include Makefile.dep diff --git a/deal.II/deal.II/Attic/examples/poisson/poisson.cc b/deal.II/deal.II/Attic/examples/poisson/poisson.cc index f4f0530b05..d6d6321fd9 100644 --- a/deal.II/deal.II/Attic/examples/poisson/poisson.cc +++ b/deal.II/deal.II/Attic/examples/poisson/poisson.cc @@ -4,13 +4,15 @@ #include #include #include -#include #include +#include +#include #include #include #include #include #include +#include #include #include @@ -20,6 +22,7 @@ extern "C" { extern TriaActiveIterator<1,CellAccessor<1> > x; extern TriaActiveIterator<2,CellAccessor<2> > y; +extern TriaRawIterator<1,DoFLineAccessor<1,LineAccessor<1> > > z; @@ -33,6 +36,16 @@ class PoissonEquation : public Equation { vector &rhs, const FEValues &fe_values, const Triangulation::cell_iterator &cell) const; + double right_hand_side (const Point &) const; +}; + + + +template +inline +double PoissonEquation::right_hand_side (const Point &) const { + // this yields as solution u=x^2+y^2+...=\vec x^2 + return -2*dim; }; @@ -45,26 +58,49 @@ void PoissonEquation<1>::assemble (dFMatrix &cell_matrix, for (unsigned int i=0; i::assemble (dFMatrix &cell_matrix, + vector &rhs, + const FEValues<2> &fe_values, + const Triangulation<2>::cell_iterator &) const { + for (unsigned int point=0; point tria; - DoFHandler<1> dof(tria); - FELinear<1> fe; - ProblemBase<1> problem(&tria, &dof); - Quadrature<1> q; +int main () { + Triangulation<2> tria; + DoFHandler<2> dof(&tria); + FELinear<2> fe; + ProblemBase<2> problem(&tria, &dof); + PoissonEquation<2> equation; + QGauss4<2> quadrature; - tria.create_hypercube(); + tria.create_unit_hypercube(); + tria.refine_global (1); dof.distribute_dofs (fe); - problem.assemble (q); + problem.assemble (equation, quadrature, fe); + + return 0; }; diff --git a/deal.II/deal.II/Todo b/deal.II/deal.II/Todo index eedc186366..42d502cd44 100644 --- a/deal.II/deal.II/Todo +++ b/deal.II/deal.II/Todo @@ -34,12 +34,19 @@ Unify CellAccessor<1> and <2> by renaming Add support for fast assemblage of mass matrices. +Make AssemblerData a local class to Assembler again if gcc2.8 supports it. DEAL: make dSMatrixStruct warn if max # of entries per row is overrun! +let dSMatrixStruct have an empty constructor to allow for member + objects being initialized later + why do loops in dFMatrix run backwards? Let an empty dfmatrix assume it has dimensions zero rather than one! + +Is dvector::operator= (double) really useful or does it make more + confusion than it helps? diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h index 305b46f1c4..a084417d2a 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.h @@ -12,6 +12,7 @@ class ostream; +class dSMatrix; class dSMatrixStruct; @@ -173,6 +174,28 @@ class ConstraintMatrix { * end of the function. */ void condense (dSMatrixStruct &sparsity) const; + + /** + * Condense a given matrix. The associated + * matrix struct should be condensed and + * compressed. + * + * The constraint matrix object must be + * closed to call this function. + */ + void condense (const dSMatrix &uncondensed, + dSMatrix &condensed) const; + + + /** + * This function does much the same as + * the above one, except that it condenses + * the matrix 'in-place'. See the general + * documentation of this class for more + * detailed information. + */ + void condense (dSMatrix &matrix) const; + /** * Print the constraint lines. Mainly for diff --git a/deal.II/deal.II/include/dofs/dof_handler.h b/deal.II/deal.II/include/dofs/dof_handler.h index 8a72bd2415..78a868da31 100644 --- a/deal.II/deal.II/include/dofs/dof_handler.h +++ b/deal.II/deal.II/include/dofs/dof_handler.h @@ -469,9 +469,10 @@ class DoFHandler : public DoFDimensionInfo { * one dimension, the functions returns * immediately after clearing the * constraint matrix. - * * For more than one dimension, the matrix - * is cleared before usage. + * is cleared before usage. The constraint + * matrix is closed anyway, no matter of the + * dimension. * * To condense a given sparsity pattern, * use #ConstraintMatrix::condense#. diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index ca758eb7e6..d23143ef7e 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -354,6 +354,10 @@ class FiniteElementBase { * Exception */ DeclException0 (ExcPureFunctionCalled); + /** + * Exception + */ + DeclException0 (ExcNotImplemented); protected: /** diff --git a/deal.II/deal.II/include/numerics/assembler.h b/deal.II/deal.II/include/numerics/assembler.h index 4f6c6ac987..a7d2e036c1 100644 --- a/deal.II/deal.II/include/numerics/assembler.h +++ b/deal.II/deal.II/include/numerics/assembler.h @@ -89,6 +89,62 @@ class Equation { +/** + * Structure to be passed upon + * construction of an assembler object + * through the iterator object. See + * \Ref{TriaRawIterator} for a discussion + * of this mechanism. + */ +template +struct AssemblerData { + /** + * Constructor. + */ + AssemblerData (const DoFHandler &dof, + const bool assemble_matrix, + const bool assemble_rhs, + const unsigned int n_rhs, + ProblemBase &problem, + const Quadrature &quadrature, + const FiniteElement &fe); + + /** + * Pointer to the dof handler object + * to be used to iterate on. + */ + const DoFHandler &dof; + /** + * Flags whether to assemble the matrix + * and right hand sides. + */ + const bool assemble_matrix, assemble_rhs; + /** + * How many right hand sides are there. + */ + const unsigned int n_rhs; + /** + * Pointer to the #ProblemBase# object + * of which the global matrices and + * vectors are to be assembled. + */ + ProblemBase &problem; + /** + * Pointer to a quadrature object to be + * used for this assemblage process. + */ + const Quadrature &quadrature; + /** + * Use this FE type for the assemblage + * process. The FE object must match that + * used to construct the #DoFHandler# + * object. + */ + const FiniteElement &fe; +}; + + + /** An #Assembler# is a specialized version of a #DoFCellAccessor# which adds @@ -97,38 +153,6 @@ class Equation { template class Assembler : public DoFCellAccessor { public: - /** - * Structure to be passed upon - * construction of an assembler object. - */ - struct AssemblerData { - /** - * Pointer to the dof handler object - * to be used to iterate on. - */ - DoFHandler *dof; - /** - * Flags whether to assemble the matrix - * and right hand sides. - */ - bool assemble_matrix, assemble_rhs; - /** - * How many right hand sides are there. - */ - unsigned int n_rhs; - /** - * Pointer to the #ProblemBase# object - * of which the global matrices and - * vectors are to be assembled. - */ - ProblemBase *problem; - /** - * Pointer to a quadrature object to be - * used for this assemblage process. - */ - Quadrature *quadrature; - }; - /** * Default constructor, unused thus not * implemented. @@ -196,7 +220,23 @@ class Assembler : public DoFCellAccessor { * of which we are to use the matrices * and vectors. */ - ProblemBase *problem; + ProblemBase &problem; + + /** + * Pointer to the finite element used for + * the assemblage process. We store a + * pointer to the finite element, not a + * copy, but this is not a problem, since + * #Assembler# objects are usually created + * and deleted by the + * #ProblemBase::assemble# routine, which + * itself gets passed the finite element; + * the lifetime of the #Assembler# object + * is thus less than the lifetime of the + * finite element object and storing + * a pointer is not risky. + */ + const FiniteElement &fe; /** * The finite element evaluated at the diff --git a/deal.II/deal.II/include/numerics/base.h b/deal.II/deal.II/include/numerics/base.h index a55cb90989..ec9f7e6548 100644 --- a/deal.II/deal.II/include/numerics/base.h +++ b/deal.II/deal.II/include/numerics/base.h @@ -6,7 +6,7 @@ #include - +#include // forward declaration template class Triangulation; @@ -17,24 +17,103 @@ template class Quadrature; +/** + Base class for user problems. This class stores the system matrix and right + hand side vectors as well as a solution vector. It initiates the assemblage + process of matrix and vectors and so on. + + This class is not extremely versatile as could certainly be. For example + it presently only supports sparse matrices and has no multigrid features. + However, all these things depend strongly on the problem and it seems + best to implement many of these things yourself. Thus, this class is more + a display of concept haw to work with deal.II. + + {\bf Assemblage} + + The #assemble# member function does the assemblage of the system matrix and + the given number of right hand sides. It does the following steps: + \begin{itemize} + \item Create sparsity pattern of the system matrix and condense it with + the constraints induced by hanging nodes. + \item Initialize an assembler object. + \item Loop over all cells and assemble matrix and vectors using the given + quadrature formula and the equation object which contains the weak + formulation of the equation. + \item Condense the system matrix with the constraints induced by hanging + nodes. + \end{itemize} + */ template class ProblemBase { public: + /** + * Constructor. Use this triangulation and + * degree of freedom object during the + * lifetime of this object. The dof + * object must refer to the given + * triangulation. The number of right hand + * sides being assembled defaults to one. + */ ProblemBase (Triangulation *tria, - DoFHandler *dof_handler); - - virtual void assemble (const Equation &equation, - const Quadrature &q); - + DoFHandler *dof_handler, + const unsigned int n_rhs=1); + + /** + * Initiate the process of assemblage of + * vectors and system matrix. Use the + * given equation object and the given + * quadrature formula. + * + * For what exactly happens here, refer to + * the general doc of this class. + */ + virtual void assemble (const Equation &equation, + const Quadrature &q, + const FiniteElement &fe); + + /** + * Exception + */ + DeclException0 (ExcDofAndTriaDontMatch); + /** + * Exception + */ + DeclException0 (ExcNoMemory); + protected: + /** + * Pointer to the triangulation to work on. + */ Triangulation *tria; + + /** + * Pointer to the degree of freedom handler + * to work with. + */ DoFHandler *dof_handler; - dSMatrix system_matrix; + /** + * Sparsity pattern of the system matrix. + */ dSMatrixStruct system_sparsity; + /** + * System matrix. + */ + dSMatrix system_matrix; + + /** + * Pointers to the right hand sides to the + * problem. Usually, one only needs one + * right hand side, but simetimes one + * wants more. + */ vector right_hand_sides; + + /** + * Solution vector. + */ dVector solution; friend class Assembler; diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index 115730e9d4..525def890f 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -235,6 +235,11 @@ void ConstraintMatrix::condense (dSMatrixStruct &sparsity) const { +void ConstraintMatrix::condense (const dSMatrix &uncondensed, + dSMatrix &condensed) const {}; +void ConstraintMatrix::condense (dSMatrix &uncondensed) const {}; + + unsigned int ConstraintMatrix::n_constraints () const { return lines.size(); diff --git a/deal.II/deal.II/source/dofs/dof_handler.cc b/deal.II/deal.II/source/dofs/dof_handler.cc index ef6bf77ba5..e88aa3bdd7 100644 --- a/deal.II/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/deal.II/source/dofs/dof_handler.cc @@ -834,6 +834,7 @@ void DoFHandler<2>::do_renumbering (const vector &new_numbers) { void DoFHandler<1>::make_constraint_matrix (ConstraintMatrix &cm) const { cm.clear (); + cm.close (); }; diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 6d70b92e3b..c2e1434359 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -26,7 +26,8 @@ FEValues::FEValues (const FiniteElement &fe, JxW_values(quadrature.n_quadrature_points, 0), quadrature_points(quadrature.n_quadrature_points, Point()), unit_quadrature_points(quadrature.n_quadrature_points, Point()), - jacobi_matrices (quadrature.n_quadrature_points) + jacobi_matrices (quadrature.n_quadrature_points, + dFMatrix(dim,dim)) { for (unsigned int i=0; i::reinit (const Triangulation<1>::cell_iterator &cell, unit_quadrature_points, jacobi_matrices, quadrature_points); + // compute gradients on real element for (unsigned int i=0; i::FELinear () : double FELinear<1>::shape_value(const unsigned int i, - const Point<1>& p) const + const Point<1>& p) const { Assert((i::shape_grad (const unsigned int i, // this function may be generalised to three or more dimensions with gcc2.8 -// you will have to change the number of vertices and the setting up of the -// jacobi matrices +// you will have to change the number of vertices void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell, const vector > &unit_points, vector &jacobians, @@ -332,6 +331,40 @@ FEQuadratic<2>::FEQuadratic () : +template +double +FEQuadratic::shape_value (const unsigned int i, + const Point &) const +{ + Assert (i +Point +FEQuadratic::shape_grad (const unsigned int i, + const Point &) const +{ + Assert (i (); +}; + + + +void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &, + const vector > &, + vector &, + vector > &) const { + Assert (false, ExcNotImplemented()); +}; + + + + FECubic<1>::FECubic () : @@ -353,3 +386,51 @@ void FECubic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell, FECubic<2>::FECubic () : FiniteElement<2> (1, 2, 4) {}; + + + +template +double +FECubic::shape_value (const unsigned int i, + const Point &) const +{ + Assert (i +Point +FECubic::shape_grad (const unsigned int i, + const Point &) const +{ + Assert (i (); +}; + + + +void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &, + const vector > &, + vector &, + vector > &) const { + Assert (false, ExcNotImplemented()); +}; + + + + + + +// explicite instantiations + +template class FELinear<1>; +template class FELinear<2>; +template class FEQuadratic<1>; +template class FEQuadratic<2>; +template class FECubic<1>; +template class FECubic<2>; + diff --git a/deal.II/deal.II/source/numerics/assembler.cc b/deal.II/deal.II/source/numerics/assembler.cc index 1a1da666a9..3b99a37c69 100644 --- a/deal.II/deal.II/source/numerics/assembler.cc +++ b/deal.II/deal.II/source/numerics/assembler.cc @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -67,18 +68,20 @@ void Equation<2>::assemble (vector &, -/* + template -Assembler::AssemblerData::AssemblerData (DoFHandler *dof, - const bool assemble_matrix, - const bool assemble_rhs, - const unsigned int n_rhs, - ProblemBase *problem, - const Quadrature *quadrature) : +AssemblerData::AssemblerData (const DoFHandler &dof, + const bool assemble_matrix, + const bool assemble_rhs, + const unsigned int n_rhs, + ProblemBase &problem, + const Quadrature &quadrature, + const FiniteElement &fe) : dof(dof), assemble_matrix(assemble_matrix), assemble_rhs(assemble_rhs), n_rhs(n_rhs), - problem(problem), quadrature(quadrature) {}; - */ + problem(problem), quadrature(quadrature), + fe(fe) {}; + @@ -88,31 +91,31 @@ Assembler::Assembler (Triangulation *tria, const int index, const void *local_data) : DoFCellAccessor (tria,level,index, - ((AssemblerData*)local_data)->dof), + &((AssemblerData*)local_data)->dof), cell_matrix (dof_handler->get_selected_fe().total_dofs), - cell_vectors (((AssemblerData*)local_data)->n_rhs, + cell_vectors (((AssemblerData*)local_data)->n_rhs, dVector(dof_handler->get_selected_fe().total_dofs)), - assemble_matrix (((AssemblerData*)local_data)->assemble_matrix), - assemble_rhs (((AssemblerData*)local_data)->assemble_rhs), - problem (((AssemblerData*)local_data)->problem), - fe_values (dof_handler->get_selected_fe(), - *((AssemblerData*)local_data)->quadrature) + assemble_matrix (((AssemblerData*)local_data)->assemble_matrix), + assemble_rhs (((AssemblerData*)local_data)->assemble_rhs), + problem (((AssemblerData*)local_data)->problem), + fe(((AssemblerData*)local_data)->fe), + fe_values (((AssemblerData*)local_data)->fe, + ((AssemblerData*)local_data)->quadrature) { - Assert (((AssemblerData*)local_data)->dof != 0, + Assert (((AssemblerData*)local_data)->n_rhs != 0, ExcInvalidData()); - Assert (((AssemblerData*)local_data)->n_rhs != 0, + Assert ((unsigned int)problem.system_matrix.m() == dof_handler->n_dofs(), ExcInvalidData()); - Assert (((AssemblerData*)local_data)->problem != 0, + Assert ((unsigned int)problem.system_matrix.n() == dof_handler->n_dofs(), ExcInvalidData()); - - Assert ((unsigned int)problem->system_matrix.m() == dof_handler->n_dofs(), + Assert (problem.right_hand_sides.size() == cell_vectors.size(), ExcInvalidData()); - Assert ((unsigned int)problem->system_matrix.n() == dof_handler->n_dofs(), - ExcInvalidData()); - Assert (problem->right_hand_sides.size() == cell_vectors.size(), + Assert (((AssemblerData*)local_data)->fe == dof_handler->get_selected_fe(), ExcInvalidData()); + + for (unsigned int i=0; iright_hand_sides[i]->n() == dof_handler->n_dofs(), + Assert ((unsigned int)problem.right_hand_sides[i]->n() == dof_handler->n_dofs(), ExcInvalidData()); }; @@ -124,7 +127,7 @@ void Assembler::assemble (const Equation &equation) { fe_values.reinit (Triangulation::cell_iterator (tria, present_level, present_index), - dof_handler->get_selected_fe()); + fe); if (assemble_matrix) // clear cell matrix @@ -170,13 +173,13 @@ void Assembler::assemble (const Equation &equation) { // distribute cell matrix for (unsigned int i=0; iget_selected_fe().total_dofs; ++i) for (unsigned int j=0; jget_selected_fe().total_dofs; ++j) - problem->system_matrix.set(dofs[i], dofs[j], cell_matrix(i,j)); + problem.system_matrix.add(dofs[i], dofs[j], cell_matrix(i,j)); if (assemble_rhs) // distribute cell vector for (unsigned int vec=0; vecget_selected_fe().total_dofs; ++j) - (*problem->right_hand_sides[vec])(dofs[j]) = cell_vectors[vec](j); + (*problem.right_hand_sides[vec])(dofs[j]) += cell_vectors[vec](j); }; @@ -189,3 +192,15 @@ template class Equation<2>; template class Assembler<1>; template class Assembler<2>; + + +template class AssemblerData<1>; +template class AssemblerData<2>; + + +template class TriaRawIterator<1,Assembler<1> >; +template class TriaIterator<1,Assembler<1> >; +template class TriaActiveIterator<1,Assembler<1> >; +template class TriaRawIterator<2,Assembler<2> >; +template class TriaIterator<2,Assembler<2> >; +template class TriaActiveIterator<2,Assembler<2> >; diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index a250604a3a..1bd73903bf 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -18,7 +18,7 @@ ProblemBase::ProblemBase (Triangulation *tria, tria(tria), dof_handler(dof), system_sparsity(1,1,1), // dummy initialisation, is later reinit'd - system_matrix(system_sparsity),// dummy initialisation, is later reinit'd + system_matrix(), // dummy initialisation, is later reinit'd right_hand_sides (n_rhs, (dVector*)0) { Assert (tria == &dof->get_tria(), ExcDofAndTriaDontMatch()); @@ -35,7 +35,8 @@ ProblemBase::ProblemBase (Triangulation *tria, template void ProblemBase::assemble (const Equation &equation, - const Quadrature &quadrature) { + const Quadrature &quadrature, + const FiniteElement &fe) { system_sparsity.reinit (dof_handler->n_dofs(), dof_handler->n_dofs(), dof_handler->max_couplings_between_dofs()); @@ -55,11 +56,12 @@ void ProblemBase::assemble (const Equation &equation, system_matrix.reinit (system_sparsity); // create assembler - Assembler::AssemblerData data (dof_handler, - true, true, //assemble matrix and rhs - right_hand_sides.size(), //one rhs - this, - &quadrature); + AssemblerData data (*dof_handler, + true, true, //assemble matrix and rhs + right_hand_sides.size(), //one rhs + *this, + quadrature, + fe); TriaActiveIterator > assembler (tria, tria->begin_active()->level(), tria->begin_active()->index(), diff --git a/tests/big-tests/Makefile b/tests/big-tests/Makefile index 719e01c5d0..e874cc83e6 100644 --- a/tests/big-tests/Makefile +++ b/tests/big-tests/Makefile @@ -5,8 +5,9 @@ INCPATH = -I../include -I../../lac/include -I../../base/include CXXFLAGS.g= -DDEBUG -g -Wall $(INCPATH) CXXFLAGS = -O3 $(INCPATH) -LIBS.g = -L../lib -lbasic.g -lgrid.g -lfe.g -lnumerics.g -llac.g -LIBS = -L../lib -lbasic -lgrid -lfe -lnumerics +LIBPATH = -L../lib -L../../lac/lib +LIBS.g = $(LIBPATH) -lbasic.g -lgrid.g -lfe.g -lnumerics.g -llac.g +LIBS = $(LIBPATH) -lbasic -lgrid -lfe -lnumerics ifeq ($(shell uname),SunOS) CXXFLAGS := -V2.7.2.3 $(CXXFLAGS) @@ -14,17 +15,22 @@ CXXFLAGS.g:= -V2.7.2.3 $(CXXFLAGS.g) endif -cc-files = grid/grid_test.cc +cc-files = grid/grid_test.cc dof/dof_test.cc poisson/poisson.cc o-files = $(cc-files:.cc=.o) h-files = $(wildcard ../include/*.h) -examples: grid/grid_test dof/dof_test +examples: $(cc-files:.cc=) +%.o : %.cc + $(CXX) $(CXXFLAGS) -o $@ -c $< -grid/grid_test: grid/grid_test.o ../lib/libgrid.a +grid/grid_test: grid/grid_test.o $(CXX) $(CXXFLAGS) -o $@ $< $(LIBS) -dof/dof_test: dof/dof_test.o ../lib/libgrid.g.a ../lib/libdeal.g.a +dof/dof_test: dof/dof_test.o + $(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS.g) -lg++ + +poisson/poisson: poisson/poisson.o $(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS.g) -lg++ @@ -52,6 +58,6 @@ clean: #to remake this file upon inclusion at the bottom #of this file. Makefile.dep: $(cc-files) $(h-files) - $(CXX) -I../include -M $(cc-files) > Makefile.dep + $(CXX) $(INCPATH) -M $(cc-files) > Makefile.dep include Makefile.dep diff --git a/tests/big-tests/Makefile.dep b/tests/big-tests/Makefile.dep index 1ce93aa62f..7dafc01152 100644 --- a/tests/big-tests/Makefile.dep +++ b/tests/big-tests/Makefile.dep @@ -12,7 +12,7 @@ grid_test.o: grid/grid_test.cc ../include/grid/tria_iterator.h \ /usr/local/lib/g++-include/new.h /usr/local/lib/g++-include/std/new.h \ /usr/local/lib/g++-include/std/cstddef.h \ /usr/local/include/g++/stl_raw_storage_iter.h \ - ../include/basic/exceptions.h ../include/grid/point.h \ + ../../base/include/base/exceptions.h ../include/grid/point.h \ ../include/grid/tria_accessor.h ../include/grid/tria.h \ /usr/local/include/g++/vector.h /usr/local/include/g++/algobase.h \ /usr/local/include/g++/pair.h /usr/local/include/g++/stl_pair.h \ @@ -31,3 +31,85 @@ grid_test.o: grid/grid_test.cc ../include/grid/tria_iterator.h \ /usr/include/math.h /usr/include/floatingpoint.h /usr/include/stdio.h \ /usr/include/sys/va_list.h /usr/include/sys/ieeefp.h \ /usr/local/include/g++/fstream.h +dof_test.o: dof/dof_test.cc ../include/grid/dof.h \ + /usr/local/include/g++/vector.h /usr/local/include/g++/algobase.h \ + /usr/local/include/g++/pair.h /usr/local/include/g++/stl_config.h \ + /usr/local/include/g++/stl_relops.h /usr/local/include/g++/stl_pair.h \ + /usr/local/include/g++/iterator.h /usr/local/include/g++/function.h \ + /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/stddef.h \ + /usr/local/include/g++/stl_function.h \ + /usr/local/include/g++/iostream.h /usr/local/include/g++/streambuf.h \ + /usr/local/include/g++/libio.h /usr/local/lib/g++-include/_G_config.h \ + /usr/local/include/g++/stl_iterator.h \ + /usr/local/include/g++/type_traits.h \ + /usr/local/include/g++/stl_construct.h \ + /usr/local/lib/g++-include/new.h /usr/local/lib/g++-include/std/new.h \ + /usr/local/lib/g++-include/std/cstddef.h \ + /usr/local/include/g++/stl_raw_storage_iter.h \ + /usr/local/include/g++/stl_algobase.h /usr/include/string.h \ + /usr/include/sys/feature_tests.h \ + /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/limits.h \ + /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/syslimits.h \ + /usr/include/limits.h /usr/include/sys/isa_defs.h \ + /usr/include/stdlib.h /usr/local/include/g++/stl_uninitialized.h \ + /usr/local/include/g++/alloc.h /usr/local/include/g++/stl_alloc.h \ + /usr/local/sparc-sun-solaris2.5.1/include/assert.h \ + /usr/local/include/g++/stl_vector.h ../include/fe/fe.h \ + ../../lac/include/lac/dfmatrix.h /usr/include/stdio.h \ + /usr/include/sys/va_list.h ../../lac/include/lac/dvector.h \ + ../../base/include/base/types.h ../../lac/include/lac/vectorbase.h \ + ../../base/include/base/exceptions.h ../../lac/include/lac/ivector.h \ + ../include/grid/point.h ../include/grid/tria.h \ + ../include/grid/tria_line.h ../include/grid/tria_quad.h \ + ../include/grid/tria_boundary.h ../include/fe/fe_lib.h \ + ../include/grid/tria_iterator.h ../include/grid/tria_accessor.h \ + ../../lac/include/lac/dsmatrix.h \ + /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/math.h \ + /usr/include/math.h /usr/include/floatingpoint.h \ + /usr/include/sys/ieeefp.h ../include/basic/parameter_handler.h \ + /usr/local/include/g++/map.h /usr/local/include/g++/tree.h \ + /usr/local/include/g++/stl_tree.h /usr/local/include/g++/stl_map.h \ + /usr/local/lib/g++-include/String.h \ + /usr/local/lib/g++-include/Regex.h ../include/grid/dof_constraints.h \ + /usr/local/include/g++/fstream.h /usr/local/include/g++/cmath +poisson.o: poisson/poisson.cc ../include/grid/dof.h \ + /usr/local/include/g++/vector.h /usr/local/include/g++/algobase.h \ + /usr/local/include/g++/pair.h /usr/local/include/g++/stl_config.h \ + /usr/local/include/g++/stl_relops.h /usr/local/include/g++/stl_pair.h \ + /usr/local/include/g++/iterator.h /usr/local/include/g++/function.h \ + /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/stddef.h \ + /usr/local/include/g++/stl_function.h \ + /usr/local/include/g++/iostream.h /usr/local/include/g++/streambuf.h \ + /usr/local/include/g++/libio.h /usr/local/lib/g++-include/_G_config.h \ + /usr/local/include/g++/stl_iterator.h \ + /usr/local/include/g++/type_traits.h \ + /usr/local/include/g++/stl_construct.h \ + /usr/local/lib/g++-include/new.h /usr/local/lib/g++-include/std/new.h \ + /usr/local/lib/g++-include/std/cstddef.h \ + /usr/local/include/g++/stl_raw_storage_iter.h \ + /usr/local/include/g++/stl_algobase.h /usr/include/string.h \ + /usr/include/sys/feature_tests.h \ + /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/limits.h \ + /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/syslimits.h \ + /usr/include/limits.h /usr/include/sys/isa_defs.h \ + /usr/include/stdlib.h /usr/local/include/g++/stl_uninitialized.h \ + /usr/local/include/g++/alloc.h /usr/local/include/g++/stl_alloc.h \ + /usr/local/sparc-sun-solaris2.5.1/include/assert.h \ + /usr/local/include/g++/stl_vector.h ../include/fe/fe.h \ + ../../lac/include/lac/dfmatrix.h /usr/include/stdio.h \ + /usr/include/sys/va_list.h ../../lac/include/lac/dvector.h \ + ../../base/include/base/types.h ../../lac/include/lac/vectorbase.h \ + ../../base/include/base/exceptions.h ../../lac/include/lac/ivector.h \ + ../include/grid/point.h ../include/grid/tria.h \ + ../include/grid/tria_line.h ../include/grid/tria_quad.h \ + ../include/grid/tria_boundary.h ../include/fe/fe_lib.h \ + ../include/grid/tria_accessor.h ../include/grid/dof_accessor.h \ + ../include/grid/tria_iterator.h ../../lac/include/lac/dsmatrix.h \ + /usr/local/lib/gcc-lib/sparc-sun-solaris2.5.1/2.8.0/include/math.h \ + /usr/include/math.h /usr/include/floatingpoint.h \ + /usr/include/sys/ieeefp.h ../include/grid/dof_constraints.h \ + ../include/basic/data_io.h /usr/local/lib/g++-include/String.h \ + /usr/local/lib/g++-include/Regex.h ../include/numerics/base.h \ + ../include/numerics/assembler.h ../include/fe/quadrature_lib.h \ + ../include/fe/quadrature.h /usr/local/include/g++/fstream.h \ + /usr/local/include/g++/cmath diff --git a/tests/big-tests/poisson/poisson.cc b/tests/big-tests/poisson/poisson.cc index f4f0530b05..d6d6321fd9 100644 --- a/tests/big-tests/poisson/poisson.cc +++ b/tests/big-tests/poisson/poisson.cc @@ -4,13 +4,15 @@ #include #include #include -#include #include +#include +#include #include #include #include #include #include +#include #include #include @@ -20,6 +22,7 @@ extern "C" { extern TriaActiveIterator<1,CellAccessor<1> > x; extern TriaActiveIterator<2,CellAccessor<2> > y; +extern TriaRawIterator<1,DoFLineAccessor<1,LineAccessor<1> > > z; @@ -33,6 +36,16 @@ class PoissonEquation : public Equation { vector &rhs, const FEValues &fe_values, const Triangulation::cell_iterator &cell) const; + double right_hand_side (const Point &) const; +}; + + + +template +inline +double PoissonEquation::right_hand_side (const Point &) const { + // this yields as solution u=x^2+y^2+...=\vec x^2 + return -2*dim; }; @@ -45,26 +58,49 @@ void PoissonEquation<1>::assemble (dFMatrix &cell_matrix, for (unsigned int i=0; i::assemble (dFMatrix &cell_matrix, + vector &rhs, + const FEValues<2> &fe_values, + const Triangulation<2>::cell_iterator &) const { + for (unsigned int point=0; point tria; - DoFHandler<1> dof(tria); - FELinear<1> fe; - ProblemBase<1> problem(&tria, &dof); - Quadrature<1> q; +int main () { + Triangulation<2> tria; + DoFHandler<2> dof(&tria); + FELinear<2> fe; + ProblemBase<2> problem(&tria, &dof); + PoissonEquation<2> equation; + QGauss4<2> quadrature; - tria.create_hypercube(); + tria.create_unit_hypercube(); + tria.refine_global (1); dof.distribute_dofs (fe); - problem.assemble (q); + problem.assemble (equation, quadrature, fe); + + return 0; };