finite element discretisation of one or more equations.
Equation objects need only provide functions which set up the cell
- matrices and the right hand side(s). These are then automatically inserted
+ matrices and the cell right hand side. These are then automatically inserted
into the global matrices and vectors.
*/
template <int dim>
/**
* Virtual function which assembles the
- * cell matrix and a specific (user
- * selectable) number of right hand sides
+ * cell matrix and the right hand side
* on a given cell.
*/
virtual void assemble (dFMatrix &cell_matrix,
- vector<dVector> &rhs,
+ dVector &rhs,
const FEValues<dim> &fe_values,
const Triangulation<dim>::cell_iterator &cell) const;
/**
* Virtual function which only assembles
- * the right hand side(s) on a given cell.
+ * the right hand side on a given cell.
*/
- virtual void assemble (vector<dVector> &rhs,
+ virtual void assemble (dVector &rhs,
const FEValues<dim> &fe_values,
const Triangulation<dim>::cell_iterator &cell) const;
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);
* 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
dFMatrix cell_matrix;
/**
- * Have room for a certain number of
- * right hand sides.
+ * Right hand side local to cell.
*/
- vector<dVector> cell_vectors;
+ dVector cell_vector;
/**
* Store whether to assemble the
/**
* Store whether to assemble the
- * right hand sides.
+ * right hand side.
*/
bool assemble_rhs;
\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.
+ \item Condense the system matrix and right hand side with the constraints
+ induced by hanging nodes.
\end{itemize}
*/
template <int dim>
* 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.
+ * triangulation.
*/
ProblemBase (Triangulation<dim> *tria,
- DoFHandler<dim> *dof_handler,
- const unsigned int n_rhs=1);
+ DoFHandler<dim> *dof_handler);
/**
* Initiate the process of assemblage of
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 storing the right hand side.
*/
- vector<dVector*> right_hand_sides;
+ dVector right_hand_side;
/**
* Solution vector.
void Equation<1>::assemble (dFMatrix &,
- vector<dVector> &,
+ dVector &,
const FEValues<1> &,
const Triangulation<1>::cell_iterator &) const {
Assert (false, ExcPureVirtualFunctionCalled());
void Equation<2>::assemble (dFMatrix &,
- vector<dVector> &,
+ dVector &,
const FEValues<2> &,
const Triangulation<2>::cell_iterator &) const {
Assert (false, ExcPureVirtualFunctionCalled());
-void Equation<1>::assemble (vector<dVector> &,
+void Equation<1>::assemble (dVector &,
const FEValues<1> &,
const Triangulation<1>::cell_iterator &) const {
Assert (false, ExcPureVirtualFunctionCalled());
-void Equation<2>::assemble (vector<dVector> &,
+void Equation<2>::assemble (dVector &,
const FEValues<2> &,
const Triangulation<2>::cell_iterator &) const {
Assert (false, ExcPureVirtualFunctionCalled());
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),
+ dof(dof),
+ assemble_matrix(assemble_matrix),
+ assemble_rhs(assemble_rhs),
+ problem(problem),
+ quadrature(quadrature),
fe(fe) {};
DoFCellAccessor<dim> (tria,level,index,
&((AssemblerData<dim>*)local_data)->dof),
cell_matrix (dof_handler->get_selected_fe().total_dofs),
- cell_vectors (((AssemblerData<dim>*)local_data)->n_rhs,
- dVector(dof_handler->get_selected_fe().total_dofs)),
+ cell_vector (dVector(dof_handler->get_selected_fe().total_dofs)),
assemble_matrix (((AssemblerData<dim>*)local_data)->assemble_matrix),
assemble_rhs (((AssemblerData<dim>*)local_data)->assemble_rhs),
problem (((AssemblerData<dim>*)local_data)->problem),
fe_values (((AssemblerData<dim>*)local_data)->fe,
((AssemblerData<dim>*)local_data)->quadrature)
{
- Assert (((AssemblerData<dim>*)local_data)->n_rhs != 0,
- ExcInvalidData());
Assert ((unsigned int)problem.system_matrix.m() == dof_handler->n_dofs(),
ExcInvalidData());
Assert ((unsigned int)problem.system_matrix.n() == dof_handler->n_dofs(),
ExcInvalidData());
- Assert (problem.right_hand_sides.size() == cell_vectors.size(),
- ExcInvalidData());
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(),
- ExcInvalidData());
+ Assert ((unsigned int)problem.right_hand_side.n() == dof_handler->n_dofs(),
+ ExcInvalidData());
};
present_level,
present_index),
fe);
-
+
+ // clear cell matrix
if (assemble_matrix)
- // clear 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)
cell_matrix(i,j) = 0;
+ // clear cell vector
if (assemble_rhs)
- // clear 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)
- cell_vectors[vec](j) = 0;
+ cell_vector.clear ();
// fill cell matrix and vector if required
if (assemble_matrix && assemble_rhs)
- equation.assemble (cell_matrix, cell_vectors, fe_values,
+ equation.assemble (cell_matrix, cell_vector, fe_values,
Triangulation<dim>::cell_iterator(tria,
present_level,
present_index));
present_index));
else
if (assemble_rhs)
- equation.assemble (cell_vectors, fe_values,
+ equation.assemble (cell_vector, fe_values,
Triangulation<dim>::cell_iterator(tria,
present_level,
present_index));
vector<int> dofs;
dof_indices (dofs);
+ // distribute cell matrix
if (assemble_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.add(dofs[i], dofs[j], cell_matrix(i,j));
+ // distribute cell vector
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);
+ for (unsigned int j=0; j<dof_handler->get_selected_fe().total_dofs; ++j)
+ problem.right_hand_side(dofs[j]) += cell_vector(j);
};
#include <grid/tria_iterator.h>
#include <basic/data_io.h>
-
-#include "../../../mia/vectormemory.h"
#include "../../../mia/control.h"
+#include "../../../mia/vectormemory.h"
#include "../../../mia/cg.h"
+
+
extern TriaActiveIterator<1,CellAccessor<1> > __dummy1233; // for gcc2.7
extern TriaActiveIterator<2,CellAccessor<2> > __dummy1234;
template <int dim>
ProblemBase<dim>::ProblemBase (Triangulation<dim> *tria,
- DoFHandler<dim> *dof,
- const unsigned int n_rhs) :
+ DoFHandler<dim> *dof) :
tria(tria),
dof_handler(dof),
system_sparsity(1,1,1), // dummy initialisation, is later reinit'd
- system_matrix(), // dummy initialisation, is later reinit'd
- right_hand_sides (n_rhs, (dVector*)0)
+ system_matrix() // dummy initialisation, is later reinit'd
{
Assert (tria == &dof->get_tria(), ExcDofAndTriaDontMatch());
-
- for (unsigned int i=0; i<n_rhs; ++i)
- {
- right_hand_sides[i] = new dVector;
- Assert (right_hand_sides[i] != 0, ExcNoMemory());
- };
};
system_sparsity.reinit (dof_handler->n_dofs(),
dof_handler->n_dofs(),
dof_handler->max_couplings_between_dofs());
- for (unsigned int i=0; i<right_hand_sides.size(); ++i)
- right_hand_sides[i]->reinit (dof_handler->n_dofs());
+ right_hand_side.reinit (dof_handler->n_dofs());
// make up sparsity pattern and
// compress with constraints
// create assembler
AssemblerData<dim> data (*dof_handler,
true, true, //assemble matrix and rhs
- right_hand_sides.size(),
*this,
quadrature,
fe);
// condense system matrix in-place
constraints.condense (system_matrix);
- // condense right hand sides in-place
- for (unsigned int i=0; i<right_hand_sides.size(); ++i)
- constraints.condense (*right_hand_sides[i]);
+
+ // condense right hand side in-place
+ constraints.condense (right_hand_side);
};
double tolerance = 1.e-16;
Control control1(max_iter,tolerance);
- PrimitiveVectorMemory<dVector> memory(right_hand_sides[0]->n());
+ PrimitiveVectorMemory<dVector> memory(right_hand_side.n());
CG<dSMatrix,dVector> cg(control1,memory);
// solve
- cg (system_matrix, solution, *right_hand_sides[0]);
+ cg (system_matrix, solution, right_hand_side);
// distribute solution
constraints.distribute (solution);
};