*/
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
--- /dev/null
+/* $Id$ */
+
+#include <fe/quadrature.h>
+
+
+template <int dim>
+Quadrature<dim>::Quadrature (const unsigned int n_q) :
+ n_quadrature_points(n_q),
+ quadrature_points (n_q, Point<dim>()),
+ weights (n_q, 0) {};
+
+
+
+template <int dim>
+const Point<dim> & Quadrature<dim>::quad_point (const unsigned int i) const {
+ Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
+ return quadrature_points[i];
+};
+
+
+
+template <int dim>
+double Quadrature<dim>::weight (const unsigned int i) const {
+ Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
+ return weights[i];
+};
+
+
+
+
+
+// explicite instantiations
+template class Quadrature<1>;
+template class Quadrature<2>;
for (unsigned int i=0; i<n_quadrature_points; ++i)
{
- quadrature_points.push_back (Point<1> (xpts[i]));
- weights.push_back (wts[i]);
+ quadrature_points[i] = Point<1>(xpts[i]);
+ weights[i] = wts[i];
};
};
for (unsigned int i=0; i<n_quadrature_points; ++i)
{
- quadrature_points.push_back (Point<1> (xpts[i]));
- weights.push_back (wts[i]);
+ quadrature_points[i] = Point<1>(xpts[i]);
+ weights[i] = wts[i];
};
};
for (unsigned int i=0; i<n_quadrature_points; ++i)
{
- quadrature_points.push_back (Point<1> (xpts[i]));
- weights.push_back (wts[i]);
+ quadrature_points[i] = Point<1>(xpts[i]);
+ weights[i] = wts[i];
};
};
for (unsigned int i=0; i<n_quadrature_points; ++i)
{
- quadrature_points.push_back (Point<1> (xpts[i]));
- weights.push_back (wts[i]);
+ quadrature_points[i] = Point<1>(xpts[i]);
+ weights[i] = wts[i];
};
};
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;
};
for (unsigned int i=0; i<n_quadrature_points; ++i)
{
- quadrature_points.push_back (Point<1> (xpts[i]));
- weights.push_back (wts[i]);
+ quadrature_points[i] = Point<1>(xpts[i]);
+ weights[i] = wts[i];
};
};
for (unsigned int i=0; i<n_quadrature_points; ++i)
{
- quadrature_points.push_back (Point<1> (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<n_quadrature_points; ++i)
+ {
+ quadrature_points[i] = Point<2>(xpts[i], ypts[i]);
+ weights[i] = wts[i];
+ };
+
+};
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)
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++
#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
#include <grid/dof.h>
#include <grid/tria.h>
#include <fe/fe_lib.h>
-#include <grid/tria_iterator.h>
#include <grid/tria_accessor.h>
+#include <grid/dof_accessor.h>
+#include <grid/tria_iterator.h>
#include <lac/dsmatrix.h>
#include <grid/dof_constraints.h>
#include <basic/data_io.h>
#include <numerics/base.h>
#include <numerics/assembler.h>
+#include <fe/quadrature_lib.h>
#include <fstream.h>
#include <cmath>
extern TriaActiveIterator<1,CellAccessor<1> > x;
extern TriaActiveIterator<2,CellAccessor<2> > y;
+extern TriaRawIterator<1,DoFLineAccessor<1,LineAccessor<1> > > z;
vector<dVector> &rhs,
const FEValues<dim> &fe_values,
const Triangulation<dim>::cell_iterator &cell) const;
+ double right_hand_side (const Point<dim> &) const;
+};
+
+
+
+template <int dim>
+inline
+double PoissonEquation<dim>::right_hand_side (const Point<dim> &) const {
+ // this yields as solution u=x^2+y^2+...=\vec x^2
+ return -2*dim;
};
for (unsigned int i=0; i<fe_values.total_dofs; ++i)
{
for (unsigned int j=0; j<fe_values.total_dofs; ++j)
- cell_matrix(i,j) = fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point) *
- fe_values.JxW(point);
- rhs[0] = fe_values.shape_value(i,point) *
- right_hand_side(fe_values.quadrature_point(point)) *
- fe_values.JxW(point);
+ cell_matrix(i,j) += fe_values.shape_grad(i,point) *
+ fe_values.shape_grad(j,point) *
+ fe_values.JxW(point);
+ rhs[0](i) += fe_values.shape_value(i,point) *
+ right_hand_side(fe_values.quadrature_point(point)) *
+ fe_values.JxW(point);
+ };
+};
+
+
+
+void PoissonEquation<2>::assemble (dFMatrix &cell_matrix,
+ vector<dVector> &rhs,
+ const FEValues<2> &fe_values,
+ const Triangulation<2>::cell_iterator &) const {
+ for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+ for (unsigned int i=0; i<fe_values.total_dofs; ++i)
+ {
+ for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+ cell_matrix(i,j) += fe_values.shape_grad(i,point) *
+ fe_values.shape_grad(j,point) *
+ fe_values.JxW(point);
+ rhs[0](i) += fe_values.shape_value(i,point) *
+ right_hand_side(fe_values.quadrature_point(point)) *
+ fe_values.JxW(point);
};
};
-void main () {
- Triangulation<1> 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;
};
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?
class ostream;
+class dSMatrix;
class dSMatrixStruct;
* 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
* 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#.
* Exception
*/
DeclException0 (ExcPureFunctionCalled);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcNotImplemented);
protected:
/**
+/**
+ * Structure to be passed upon
+ * construction of an assembler object
+ * through the iterator object. See
+ * \Ref{TriaRawIterator} for a discussion
+ * of this mechanism.
+ */
+template <int dim>
+struct AssemblerData {
+ /**
+ * Constructor.
+ */
+ AssemblerData (const DoFHandler<dim> &dof,
+ const bool assemble_matrix,
+ const bool assemble_rhs,
+ const unsigned int n_rhs,
+ ProblemBase<dim> &problem,
+ const Quadrature<dim> &quadrature,
+ const FiniteElement<dim> &fe);
+
+ /**
+ * Pointer to the dof handler object
+ * to be used to iterate on.
+ */
+ const DoFHandler<dim> &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<dim> &problem;
+ /**
+ * Pointer to a quadrature object to be
+ * used for this assemblage process.
+ */
+ const Quadrature<dim> &quadrature;
+ /**
+ * Use this FE type for the assemblage
+ * process. The FE object must match that
+ * used to construct the #DoFHandler#
+ * object.
+ */
+ const FiniteElement<dim> &fe;
+};
+
+
+
/**
An #Assembler# is a specialized version of a #DoFCellAccessor# which adds
template <int dim>
class Assembler : public DoFCellAccessor<dim> {
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<dim> *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<dim> *problem;
- /**
- * Pointer to a quadrature object to be
- * used for this assemblage process.
- */
- Quadrature<dim> *quadrature;
- };
-
/**
* Default constructor, unused thus not
* implemented.
* of which we are to use the matrices
* and vectors.
*/
- ProblemBase<dim> *problem;
+ ProblemBase<dim> &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<dim> &fe;
/**
* The finite element evaluated at the
#include <lac/dsmatrix.h>
-
+#include <base/exceptions.h>
// forward declaration
template <int dim> class Triangulation;
+/**
+ 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 <int dim>
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<dim> *tria,
- DoFHandler<dim> *dof_handler);
-
- virtual void assemble (const Equation<dim> &equation,
- const Quadrature<dim> &q);
-
+ DoFHandler<dim> *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<dim> &equation,
+ const Quadrature<dim> &q,
+ const FiniteElement<dim> &fe);
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcDofAndTriaDontMatch);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcNoMemory);
+
protected:
+ /**
+ * Pointer to the triangulation to work on.
+ */
Triangulation<dim> *tria;
+
+ /**
+ * Pointer to the degree of freedom handler
+ * to work with.
+ */
DoFHandler<dim> *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<dVector*> right_hand_sides;
+
+ /**
+ * Solution vector.
+ */
dVector solution;
friend class Assembler<dim>;
+void ConstraintMatrix::condense (const dSMatrix &uncondensed,
+ dSMatrix &condensed) const {};
+void ConstraintMatrix::condense (dSMatrix &uncondensed) const {};
+
+
unsigned int ConstraintMatrix::n_constraints () const {
return lines.size();
void DoFHandler<1>::make_constraint_matrix (ConstraintMatrix &cm) const {
cm.clear ();
+ cm.close ();
};
JxW_values(quadrature.n_quadrature_points, 0),
quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
unit_quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
- jacobi_matrices (quadrature.n_quadrature_points)
+ jacobi_matrices (quadrature.n_quadrature_points,
+ dFMatrix(dim,dim))
{
for (unsigned int i=0; i<fe.total_dofs; ++i)
for (unsigned int j=0; j<n_quadrature_points; ++j)
unit_quadrature_points,
jacobi_matrices,
quadrature_points);
+
// compute gradients on real element
for (unsigned int i=0; i<fe.total_dofs; ++i)
for (unsigned int j=0; j<n_quadrature_points; ++j)
double
FELinear<1>::shape_value(const unsigned int i,
- const Point<1>& p) const
+ const Point<1>& p) const
{
Assert((i<total_dofs), ExcInvalidIndex(i));
switch (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<Point<2> > &unit_points,
vector<dFMatrix> &jacobians,
+template <int dim>
+double
+FEQuadratic<dim>::shape_value (const unsigned int i,
+ const Point<dim> &) const
+{
+ Assert (i<total_dofs, ExcInvalidIndex(i));
+ Assert (false, ExcNotImplemented());
+ return 0.;
+};
+
+
+
+template <int dim>
+Point<dim>
+FEQuadratic<dim>::shape_grad (const unsigned int i,
+ const Point<dim> &) const
+{
+ Assert (i<total_dofs, ExcInvalidIndex(i));
+ Assert (false, ExcNotImplemented());
+ return Point<dim> ();
+};
+
+
+
+void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
+ const vector<Point<2> > &,
+ vector<dFMatrix> &,
+ vector<Point<2> > &) const {
+ Assert (false, ExcNotImplemented());
+};
+
+
+
+
FECubic<1>::FECubic () :
FECubic<2>::FECubic () :
FiniteElement<2> (1, 2, 4) {};
+
+
+
+template <int dim>
+double
+FECubic<dim>::shape_value (const unsigned int i,
+ const Point<dim> &) const
+{
+ Assert (i<total_dofs, ExcInvalidIndex(i));
+ Assert (false, ExcNotImplemented());
+ return 0.;
+};
+
+
+
+template <int dim>
+Point<dim>
+FECubic<dim>::shape_grad (const unsigned int i,
+ const Point<dim> &) const
+{
+ Assert (i<total_dofs, ExcInvalidIndex(i));
+ Assert (false, ExcNotImplemented());
+ return Point<dim> ();
+};
+
+
+
+void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
+ const vector<Point<2> > &,
+ vector<dFMatrix> &,
+ vector<Point<2> > &) 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>;
+
#include <numerics/assembler.h>
#include <numerics/base.h>
#include <grid/tria_iterator.h>
+#include <grid/tria_iterator.templates.h>
#include <lac/dfmatrix.h>
#include <lac/dvector.h>
-/*
+
template <int dim>
-Assembler<dim>::AssemblerData::AssemblerData (DoFHandler<dim> *dof,
- const bool assemble_matrix,
- const bool assemble_rhs,
- const unsigned int n_rhs,
- ProblemBase<dim> *problem,
- const Quadrature<dim> *quadrature) :
+AssemblerData<dim>::AssemblerData (const DoFHandler<dim> &dof,
+ const bool assemble_matrix,
+ const bool assemble_rhs,
+ const unsigned int n_rhs,
+ ProblemBase<dim> &problem,
+ const Quadrature<dim> &quadrature,
+ const FiniteElement<dim> &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) {};
+
const int index,
const void *local_data) :
DoFCellAccessor<dim> (tria,level,index,
- ((AssemblerData*)local_data)->dof),
+ &((AssemblerData<dim>*)local_data)->dof),
cell_matrix (dof_handler->get_selected_fe().total_dofs),
- cell_vectors (((AssemblerData*)local_data)->n_rhs,
+ cell_vectors (((AssemblerData<dim>*)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<dim>*)local_data)->assemble_matrix),
+ assemble_rhs (((AssemblerData<dim>*)local_data)->assemble_rhs),
+ problem (((AssemblerData<dim>*)local_data)->problem),
+ fe(((AssemblerData<dim>*)local_data)->fe),
+ fe_values (((AssemblerData<dim>*)local_data)->fe,
+ ((AssemblerData<dim>*)local_data)->quadrature)
{
- Assert (((AssemblerData*)local_data)->dof != 0,
+ Assert (((AssemblerData<dim>*)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<dim>*)local_data)->fe == dof_handler->get_selected_fe(),
ExcInvalidData());
+
+
for (unsigned int i=0; i<cell_vectors.size(); ++i)
- Assert ((unsigned int)problem->right_hand_sides[i]->n() == dof_handler->n_dofs(),
+ Assert ((unsigned int)problem.right_hand_sides[i]->n() == dof_handler->n_dofs(),
ExcInvalidData());
};
fe_values.reinit (Triangulation<dim>::cell_iterator (tria,
present_level,
present_index),
- dof_handler->get_selected_fe());
+ fe);
if (assemble_matrix)
// clear cell matrix
// distribute cell matrix
for (unsigned int i=0; i<dof_handler->get_selected_fe().total_dofs; ++i)
for (unsigned int j=0; j<dof_handler->get_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; vec<cell_vectors.size(); ++vec)
for (unsigned int j=0; j<dof_handler->get_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);
};
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> >;
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());
template <int dim>
void ProblemBase<dim>::assemble (const Equation<dim> &equation,
- const Quadrature<dim> &quadrature) {
+ const Quadrature<dim> &quadrature,
+ const FiniteElement<dim> &fe) {
system_sparsity.reinit (dof_handler->n_dofs(),
dof_handler->n_dofs(),
dof_handler->max_couplings_between_dofs());
system_matrix.reinit (system_sparsity);
// create assembler
- Assembler<dim>::AssemblerData data (dof_handler,
- true, true, //assemble matrix and rhs
- right_hand_sides.size(), //one rhs
- this,
- &quadrature);
+ AssemblerData<dim> data (*dof_handler,
+ true, true, //assemble matrix and rhs
+ right_hand_sides.size(), //one rhs
+ *this,
+ quadrature,
+ fe);
TriaActiveIterator<dim, Assembler<dim> > assembler (tria,
tria->begin_active()->level(),
tria->begin_active()->index(),
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)
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++
#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
/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 \
/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
#include <grid/dof.h>
#include <grid/tria.h>
#include <fe/fe_lib.h>
-#include <grid/tria_iterator.h>
#include <grid/tria_accessor.h>
+#include <grid/dof_accessor.h>
+#include <grid/tria_iterator.h>
#include <lac/dsmatrix.h>
#include <grid/dof_constraints.h>
#include <basic/data_io.h>
#include <numerics/base.h>
#include <numerics/assembler.h>
+#include <fe/quadrature_lib.h>
#include <fstream.h>
#include <cmath>
extern TriaActiveIterator<1,CellAccessor<1> > x;
extern TriaActiveIterator<2,CellAccessor<2> > y;
+extern TriaRawIterator<1,DoFLineAccessor<1,LineAccessor<1> > > z;
vector<dVector> &rhs,
const FEValues<dim> &fe_values,
const Triangulation<dim>::cell_iterator &cell) const;
+ double right_hand_side (const Point<dim> &) const;
+};
+
+
+
+template <int dim>
+inline
+double PoissonEquation<dim>::right_hand_side (const Point<dim> &) const {
+ // this yields as solution u=x^2+y^2+...=\vec x^2
+ return -2*dim;
};
for (unsigned int i=0; i<fe_values.total_dofs; ++i)
{
for (unsigned int j=0; j<fe_values.total_dofs; ++j)
- cell_matrix(i,j) = fe_values.shape_grad(i,point) *
- fe_values.shape_grad(j,point) *
- fe_values.JxW(point);
- rhs[0] = fe_values.shape_value(i,point) *
- right_hand_side(fe_values.quadrature_point(point)) *
- fe_values.JxW(point);
+ cell_matrix(i,j) += fe_values.shape_grad(i,point) *
+ fe_values.shape_grad(j,point) *
+ fe_values.JxW(point);
+ rhs[0](i) += fe_values.shape_value(i,point) *
+ right_hand_side(fe_values.quadrature_point(point)) *
+ fe_values.JxW(point);
+ };
+};
+
+
+
+void PoissonEquation<2>::assemble (dFMatrix &cell_matrix,
+ vector<dVector> &rhs,
+ const FEValues<2> &fe_values,
+ const Triangulation<2>::cell_iterator &) const {
+ for (unsigned int point=0; point<fe_values.n_quadrature_points; ++point)
+ for (unsigned int i=0; i<fe_values.total_dofs; ++i)
+ {
+ for (unsigned int j=0; j<fe_values.total_dofs; ++j)
+ cell_matrix(i,j) += fe_values.shape_grad(i,point) *
+ fe_values.shape_grad(j,point) *
+ fe_values.JxW(point);
+ rhs[0](i) += fe_values.shape_value(i,point) *
+ right_hand_side(fe_values.quadrature_point(point)) *
+ fe_values.JxW(point);
};
};
-void main () {
- Triangulation<1> 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;
};