const unsigned int n_iterations = solve ();
std::cout << " Solver converged in " << n_iterations
- << " iterations." << std::endl;
+ << " iterations." << std::endl;
output_results ();
// @sect3{Generate some output for a given mesh}
template<int dim>
-void mesh_info(Triangulation<dim> & tria, const char * filename)
+void mesh_info(Triangulation<dim> &tria, const char *filename)
{
- // some general info:
+ // some general info:
std::cout << "Mesh info:" << std::endl
- << " dimension: " << dim << std::endl
- << " no. of cells: " << tria.n_active_cells() << std::endl;
+ << " dimension: " << dim << std::endl
+ << " no. of cells: " << tria.n_active_cells() << std::endl;
// loop over all the cells and find how often each boundary indicator is used:
{
std::map<unsigned int, unsigned int> boundary_count;
Triangulation<2>::active_cell_iterator
- cell = tria.begin_active(),
- endc = tria.end();
+ cell = tria.begin_active(),
+ endc = tria.end();
for (; cell!=endc; ++cell)
{
- for (unsigned int face=0; face<GeometryInfo<2>::faces_per_cell; ++face)
- {
- if (cell->face(face)->at_boundary())
- boundary_count[cell->face(face)->boundary_indicator()]++;
- }
+ for (unsigned int face=0; face<GeometryInfo<2>::faces_per_cell; ++face)
+ {
+ if (cell->face(face)->at_boundary())
+ boundary_count[cell->face(face)->boundary_indicator()]++;
+ }
}
-
+
std::cout << " boundary indicators: ";
for (std::map<unsigned int, unsigned int>::iterator it=boundary_count.begin();
- it!=boundary_count.end();
- ++it)
+ it!=boundary_count.end();
+ ++it)
{
- std::cout << it->first << "(" << it->second << ") ";
+ std::cout << it->first << "(" << it->second << ") ";
}
std::cout << std::endl;
}
-
+
// Now we want to write a graphical representation of the mesh to an output
// file.
std::ofstream out (filename);
gridin.attach_triangulation(triangulation);
std::ifstream f("untitled.msh");
gridin.read_msh(f);
-
+
mesh_info(triangulation, "grid-1.eps");
-
+
}
Triangulation<2> tria1, tria2;
-
+
GridGenerator::hyper_cube_with_cylindrical_hole (tria1, 0.25, 1.0);
std::vector< unsigned int > repetitions(2);
repetitions[0]=3;
repetitions[1]=2;
GridGenerator::subdivided_hyper_rectangle (tria2, repetitions,
- Point<2>(1.0,-1.0), Point<2>(4.0,1.0));
+ Point<2>(1.0,-1.0), Point<2>(4.0,1.0));
+
+ GridGenerator::merge_triangulations (tria1, tria2, triangulation);
- GridGenerator::merge_triangulations (tria1, tria2, triangulation);
-
mesh_info(triangulation, "grid-2.eps");
}
GridGenerator::hyper_cube_with_cylindrical_hole (triangulation, 0.25, 1.0);
Triangulation<2>::active_cell_iterator
- cell = triangulation.begin_active(),
- endc = triangulation.end();
+ cell = triangulation.begin_active(),
+ endc = triangulation.end();
for (; cell!=endc; ++cell)
{
- for (unsigned int i=0;i<GeometryInfo<2>::vertices_per_cell;++i)
- {
- Point<2> & v = cell->vertex(i);
- if (std::abs(v(1)-1.0)<1e-5)
- v(1)+=0.5;
- }
-
+ for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_cell; ++i)
+ {
+ Point<2> &v = cell->vertex(i);
+ if (std::abs(v(1)-1.0)<1e-5)
+ v(1)+=0.5;
+ }
+
}
-
+
const HyperBallBoundary<2> boundary_description(Point<2>(0,0),0.25);
triangulation.set_boundary (1, boundary_description);
triangulation.refine_global(2);
mesh_info(triangulation, "grid-3.eps");
- // remove boundary object from Triangulation again:
- triangulation.set_boundary (1);
+ // remove boundary object from Triangulation again:
+ triangulation.set_boundary (1);
}
// we need operator- for the iterator type anyway, do the
// test as follows, rather than via 'last >= first'
Assert (last - first >= 0,
- ExcMessage ("The given iterators do not satisfy the proper ordering."));
+ ExcMessage ("The given iterators do not satisfy the proper ordering."));
unsigned int len = last-first;
* constructor.
*/
void initialize_embedding ();
-
+
/**
* Initialize the restriction
* matrices. Called from the
static
void
parallelogram(Triangulation<dim> &tria,
- const Point<dim> (&corners)[dim],
- const bool colorize=false);
+ const Point<dim> (&corners)[dim],
+ const bool colorize=false);
/**
static
void
parallelogram(Triangulation<dim> &tria,
- const Tensor<2,dim> &corners,
- const bool colorize=false) DEAL_II_DEPRECATED;
+ const Tensor<2,dim> &corners,
+ const bool colorize=false) DEAL_II_DEPRECATED;
/**
* A parallelepiped. The first corner point is the origin. The
* function.
*/
template <int dim>
- static
- void
- parallelepiped (Triangulation<dim> &tria,
- const Point<dim> (&corners) [dim],
- const bool colorize = false);
-
+ static
+ void
+ parallelepiped (Triangulation<dim> &tria,
+ const Point<dim> (&corners) [dim],
+ const bool colorize = false);
+
/**
* A subdivided parallelepiped. The first corner point is the
* origin. The <tt>dim</tt> adjacent points are vectors describing
* function.
*/
template <int dim>
- static
- void
- subdivided_parallelepiped (Triangulation<dim> &tria,
- const unsigned int n_subdivisions,
- const Point<dim> (&corners) [dim],
- const bool colorize = false);
+ static
+ void
+ subdivided_parallelepiped (Triangulation<dim> &tria,
+ const unsigned int n_subdivisions,
+ const Point<dim> (&corners) [dim],
+ const bool colorize = false);
/**
* A subdivided parallelepiped, ie. the same as above, but where the
* function.
*/
template <int dim>
- static
- void
- subdivided_parallelepiped (Triangulation<dim> &tria,
- const unsigned int ( n_subdivisions) [dim],
- const Point<dim> (&corners) [dim],
- const bool colorize = false);
+ static
+ void
+ subdivided_parallelepiped (Triangulation<dim> &tria,
+ const unsigned int ( n_subdivisions) [dim],
+ const Point<dim> (&corners) [dim],
+ const bool colorize = false);
/**
--block;
return std::pair<unsigned int,unsigned int>(block,
- i-start_indices[block]);
+ i-start_indices[block]);
}
(const types::global_dof_index global_index) const
{
const std::pair<unsigned int,unsigned int> local_index
- = block_indices.global_to_local (global_index);
-
+ = block_indices.global_to_local (global_index);
+
return components[local_index.first].in_local_range (global_index);
}
// parallel vectors, this can only work if we can put a compress()
// in between, but we don't want to call compress() twice per entry
for (std::vector<ConstraintLine>::const_iterator
- constraint_line = lines.begin();
+ constraint_line = lines.begin();
constraint_line!=lines.end(); ++constraint_line)
{
// in case the constraint is
const typename VectorType::value_type old_value = vec(constraint_line->line);
for (unsigned int q=0; q!=constraint_line->entries.size(); ++q)
- if (vec.in_local_range(constraint_line->entries[q].first) == true)
+ if (vec.in_local_range(constraint_line->entries[q].first) == true)
vec(constraint_line->entries[q].first)
+= (static_cast<typename VectorType::value_type>
(old_value) *
vec.compress();
for (std::vector<ConstraintLine>::const_iterator
- constraint_line = lines.begin();
+ constraint_line = lines.begin();
constraint_line!=lines.end(); ++constraint_line)
- if (vec.in_local_range(constraint_line->line) == true)
+ if (vec.in_local_range(constraint_line->line) == true)
vec(constraint_line->line) = 0.;
}
// copy entries if column will not
// be condensed away, distribute
// otherwise
- for (typename SparseMatrix<number>::const_iterator
- p = uncondensed.begin(row);
- p != uncondensed.end(row); ++p)
+ for (typename SparseMatrix<number>::const_iterator
+ p = uncondensed.begin(row);
+ p != uncondensed.end(row); ++p)
if (new_line[p->column()] != -1)
condensed.add (new_line[row],
- new_line[p->column()],
+ new_line[p->column()],
p->value());
else
{
else
// line must be distributed
{
- for (typename SparseMatrix<number>::const_iterator
- p = uncondensed.begin(row);
- p != uncondensed.end(row); ++p)
+ for (typename SparseMatrix<number>::const_iterator
+ p = uncondensed.begin(row);
+ p != uncondensed.end(row); ++p)
// for each column: distribute
if (new_line[p->column()] != -1)
// column is not constrained
if (use_vectors == true)
for (unsigned int q=0; q!=next_constraint->entries.size(); ++q)
condensed_vector (new_line[next_constraint->entries[q].first])
- -= p->value() *
+ -= p->value() *
next_constraint->entries[q].second *
c->inhomogeneity;
}
&sparsity_struct[row_start[row]];
typename SparseMatrix<number>::iterator
val_ptr = (sparsity.row_length(row) == 0 ?
- sparse_matrix->end() :
+ sparse_matrix->end() :
sparse_matrix->begin(row));
const bool optimize_diagonal = sparsity.optimize_diagonal();
unsigned int counter = optimize_diagonal;
/**
* A function that can be used in SLEPc as a callback to check on
- * convergence.
+ * convergence.
*
* @note This function is redundant.
*/
- static
+ static
int
convergence_test (EPS eps,
- PetscScalar kr,
- PetscScalar ki,
- PetscReal residual_norm,
- PetscReal *estimated_error,
+ PetscScalar kr,
+ PetscScalar ki,
+ PetscReal residual_norm,
+ PetscReal *estimated_error,
void *solver_control);
-
+
/**
* Objects of this type are explicitly created, but are destroyed
* when the surrounding solver object goes out of scope, or when
* solver. Usage: Largest values of spectrum only, all problem
* types, complex.
*
- * @ingroup SLEPcWrappers
+ * @ingroup SLEPcWrappers
*
* @author Toby D. Young 2010
*/
* An implementation of the solver interface using the SLEPc
* Davidson solver. Usage (incomplete/untested): All problem types.
*
- * @ingroup SLEPcWrappers
+ * @ingroup SLEPcWrappers
*
* @author Toby D. Young 2010
*/
* change that.
*/
SolverGeneralizedDavidson (SolverControl &cn,
- const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
- const AdditionalData &data = AdditionalData());
+ const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
+ const AdditionalData &data = AdditionalData());
protected:
* An implementation of the solver interface using the SLEPc
* Jacobi-Davidson solver. Usage: All problem types.
*
- * @ingroup SLEPcWrappers
+ * @ingroup SLEPcWrappers
*
* @author Toby D. Young 2013
*/
* change that.
*/
SolverJacobiDavidson (SolverControl &cn,
- const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
- const AdditionalData &data = AdditionalData());
+ const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
+ const AdditionalData &data = AdditionalData());
protected:
{
unsigned int n_converged = 0;
- // Set the matrices of the problem
+ // Set the matrices of the problem
set_matrices (A);
- // and solve
+ // and solve
solve (n_eigenvectors, &n_converged);
-
+
if (n_converged > n_eigenvectors)
n_converged = n_eigenvectors;
AssertThrow (n_converged == n_eigenvectors,
{
unsigned int n_converged = 0;
- // Set the matrices of the problem
+ // Set the matrices of the problem
set_matrices (A, B);
// and solve
// and recreate
own_sparsity = new SparsityPattern(matrix_sparsity,
- matrix_sparsity.max_entries_per_row()
- +2*data.extra_off_diagonals,
- data.extra_off_diagonals);
+ matrix_sparsity.max_entries_per_row()
+ +2*data.extra_off_diagonals,
+ data.extra_off_diagonals);
own_sparsity->compress();
sparsity_pattern_to_use = own_sparsity;
}
// non-diagonal element in this row
Assert (this->cols->optimize_diagonal(), ExcNotImplemented());
typename SparseMatrix<number>::iterator
- diagonal_element = this->begin(row);
+ diagonal_element = this->begin(row);
number rowsum = 0;
for (typename SparseMatrix<number>::iterator
- p = diagonal_element + 1;
- p != this->end(row); ++p)
+ p = diagonal_element + 1;
+ p != this->end(row); ++p)
rowsum += std::fabs(p->value());
diagonal_element->value() += this->get_strengthen_diagonal (rowsum, row) *
- rowsum;
+ rowsum;
}
}
* of the current kind.
*/
typedef
- const Accessor<number,Constness> & value_type;
+ const Accessor<number,Constness> &value_type;
/**
* Constructor. Create an iterator into the matrix @p matrix for the given
* doing here.
*/
number nth_entry_in_row (const unsigned int row,
- const unsigned int index) const;
+ const unsigned int index) const;
// make all other sparse matrices friends
template <typename somenumber> friend class SparseMatrix;
inline
number
SparseMatrix<number>::nth_entry_in_row (const unsigned int row,
- const unsigned int index) const
+ const unsigned int index) const
{
Assert(row<cols->rows, ExcIndexRange(row,0,cols->rows));
Assert(index<cols->row_length(row),
Accessor<number, false>::Reference::operator number() const
{
return accessor->matrix->nth_entry_in_row(accessor->a_row,
- accessor->a_index);
+ accessor->a_index);
}
const unsigned int this_position
= (*this != (*this)->get_matrix().end()
- ?
- sparsity.get_rowstart_indices()[(*this)->row()] + (*this)->index()
- :
- sparsity.get_rowstart_indices()[sparsity.n_rows()]);
+ ?
+ sparsity.get_rowstart_indices()[(*this)->row()] + (*this)->index()
+ :
+ sparsity.get_rowstart_indices()[sparsity.n_rows()]);
const unsigned int other_position
= (other != (*this)->get_matrix().end()
- ?
- sparsity.get_rowstart_indices()[other->row()] + other->index()
- :
- sparsity.get_rowstart_indices()[sparsity.n_rows()]);
+ ?
+ sparsity.get_rowstart_indices()[other->row()] + other->index()
+ :
+ sparsity.get_rowstart_indices()[sparsity.n_rows()]);
return (this_position - other_position);
}
// work on the lower left part of the matrix. we know
// it's symmetric, so we can work with this alone
for (typename SparseMatrix<somenumber>::const_iterator
- p = matrix.begin(row)+1;
- (p != matrix.end(row)) && (p->column() < p->row());
- ++p)
+ p = matrix.begin(row)+1;
+ (p != matrix.end(row)) && (p->column() < p->row());
+ ++p)
temp1 += p->value() / diag[p->column()] * inner_sums[p->column()];
Assert(temp-temp1 > 0, ExcStrengthenDiagonalTooSmall());
const unsigned int *const first_after_diagonal = this->prebuilt_lower_bound[row];
number rowsum = 0;
for (typename SparseMatrix<number>::const_iterator
- p = this->begin(row)+1;
- p != this->end(row); ++p)
+ p = this->begin(row)+1;
+ p != this->end(row); ++p)
if (p->column() > p->row())
rowsum += p->value();
// get start of this row. skip
// the diagonal element
for (typename SparseMatrix<number>::const_iterator
- p = this->begin(row)+1;
- (p != this->end(row)) && (p->column() < row);
- ++p)
+ p = this->begin(row)+1;
+ (p != this->end(row)) && (p->column() < row);
+ ++p)
dst(row) -= p->value() * dst(p->column());
dst(row) *= inv_diag[row];
{
// get end of this row
for (typename SparseMatrix<number>::const_iterator
- p = this->begin(row)+1;
- p != this->end(row);
- ++p)
+ p = this->begin(row)+1;
+ p != this->end(row);
+ ++p)
if (p->column() > static_cast<unsigned int>(row))
dst(row) -= p->value() * dst(p->column());
++cell)
{
const unsigned int
- *it = begin_indices (renumbering[cell]),
- *end_cell = end_indices (renumbering[cell]);
+ *it = begin_indices (renumbering[cell]),
+ *end_cell = end_indices (renumbering[cell]);
for ( ; it != end_cell; ++it)
if (row_lengths[*it]>0)
connectivity_dof.add(*it, block);
else
{
const unsigned int
- *it = begin_indices (block),
- *end_cell = end_indices (block);
+ *it = begin_indices (block),
+ *end_cell = end_indices (block);
for ( ; it != end_cell; ++it)
if (row_lengths[*it]>0)
connectivity_dof.add(*it, block);
{
// apply renumbering when we do blocking
const unsigned int
- *it = begin_indices (renumbering[cell]),
- *end_cell = end_indices (renumbering[cell]);
+ *it = begin_indices (renumbering[cell]),
+ *end_cell = end_indices (renumbering[cell]);
for ( ; it != end_cell; ++it)
if (row_lengths[*it] > 0)
{
FPArrayComparator<Number> comparator(jacobian_size);
typedef Tensor<1,VectorizedArray<Number>::n_array_elements,Number> VEC_ARRAY;
std::map<Tensor<1,dim,VEC_ARRAY>, unsigned int,
- FPArrayComparator<Number> > cartesians(comparator);
+ FPArrayComparator<Number> > cartesians(comparator);
std::map<Tensor<2,dim,VEC_ARRAY>, unsigned int,
- FPArrayComparator<Number> > affines(comparator);
+ FPArrayComparator<Number> > affines(comparator);
// loop over all cells
for (unsigned int cell=0; cell<n_macro_cells; ++cell)
{
cartesian_data.resize(cartesians.size());
for (typename std::map<Tensor<1,dim,VEC_ARRAY>,
- unsigned int, FPArrayComparator<Number> >::iterator
- it = cartesians.begin(); it != cartesians.end(); ++it)
+ unsigned int, FPArrayComparator<Number> >::iterator
+ it = cartesians.begin(); it != cartesians.end(); ++it)
{
VectorizedArray<Number> det = make_vectorized_array<Number>(1.);
for (unsigned int d=0; d<dim; ++d)
}
affine_data.resize(affines.size());
for (typename std::map<Tensor<2,dim,VEC_ARRAY>,
- unsigned int, FPArrayComparator<Number> >::iterator
- it = affines.begin(); it != affines.end(); ++it)
+ unsigned int, FPArrayComparator<Number> >::iterator
+ it = affines.begin(); it != affines.end(); ++it)
{
Tensor<2,dim,VectorizedArray<Number> > jac;
for (unsigned int d=0; d<dim; ++d)
}
}
} // end loop over all entries in vectorization (vectorization_length
- // cells)
+ // cells)
// set information for next cell
for (unsigned int j=0; j<vectorization_length; ++j)
{
AssertIndexRange (row, constraint_pool_row_index.size()-1);
return constraint_pool_data.empty() ? 0 :
- &constraint_pool_data[0] + constraint_pool_row_index[row];
+ &constraint_pool_data[0] + constraint_pool_row_index[row];
}
{
AssertIndexRange (row, constraint_pool_row_index.size()-1);
return constraint_pool_data.empty() ? 0 :
- &constraint_pool_data[0] + constraint_pool_row_index[row+1];
+ &constraint_pool_data[0] + constraint_pool_row_index[row+1];
}
// set up sparsity structure
CompressedSparsityPattern c_sparsity(dof.n_boundary_dofs (boundary_functions),
- dof.n_boundary_dofs (boundary_functions));
+ dof.n_boundary_dofs (boundary_functions));
DoFTools::make_boundary_sparsity_pattern (dof,
boundary_functions,
dof_to_boundary_mapping,
{
int n_jobs=1;
(void) MPI_Comm_size (mpi_communicator, &n_jobs);
-
+
return n_jobs;
}
{
int rank=0;
(void) MPI_Comm_rank (mpi_communicator, &rank);
-
+
return rank;
}
-
-
+
+
MPI_Comm duplicate_communicator (const MPI_Comm &mpi_communicator)
{
MPI_Comm new_communicator;
MPI_Comm_dup (mpi_communicator, &new_communicator);
return new_communicator;
}
-
-
+
+
std::vector<unsigned int>
compute_point_to_point_communication_pattern (const MPI_Comm &mpi_comm,
const std::vector<unsigned int> &destinations)
{
unsigned int myid = Utilities::MPI::this_mpi_process(mpi_comm);
unsigned int n_procs = Utilities::MPI::n_mpi_processes(mpi_comm);
-
+
for (unsigned int i=0; i<destinations.size(); ++i)
{
Assert (destinations[i] < n_procs,
Assert (destinations[i] != myid,
ExcMessage ("There is no point in communicating with ourselves."));
}
-
-
+
+
// let all processors
// communicate the maximal
// number of destinations they
// have
const unsigned int max_n_destinations
= Utilities::MPI::max (destinations.size(), mpi_comm);
-
+
// now that we know the number
// of data packets every
// processor wants to send, set
numbers::invalid_unsigned_int);
std::copy (destinations.begin(), destinations.end(),
my_destinations.begin());
-
+
// now exchange these (we could
// communicate less data if we
// used MPI_Allgatherv, but
MPI_Allgather (&my_destinations[0], max_n_destinations, MPI_UNSIGNED,
&all_destinations[0], max_n_destinations, MPI_UNSIGNED,
mpi_comm);
-
+
// now we know who is going to
// communicate with
// whom. collect who is going
else if (all_destinations[i*max_n_destinations + j] ==
numbers::invalid_unsigned_int)
break;
-
+
return origins;
}
-
-
+
+
namespace
{
// custom MIP_Op for
{
const MinMaxAvg *in_lhs = static_cast<const MinMaxAvg *>(in_lhs_);
MinMaxAvg *inout_rhs = static_cast<MinMaxAvg *>(inout_rhs_);
-
+
Assert(*len==1, ExcInternalError());
-
+
inout_rhs->sum += in_lhs->sum;
if (inout_rhs->min>in_lhs->min)
{
if (inout_rhs->min_index > in_lhs->min_index)
inout_rhs->min_index = in_lhs->min_index;
}
-
+
if (inout_rhs->max < in_lhs->max)
{
inout_rhs->max = in_lhs->max;
}
}
}
-
-
-
+
+
+
MinMaxAvg
min_max_avg(const double my_value,
const MPI_Comm &mpi_communicator)
{
MinMaxAvg result;
-
+
const unsigned int my_id
= dealii::Utilities::MPI::this_mpi_process(mpi_communicator);
const unsigned int numproc
= dealii::Utilities::MPI::n_mpi_processes(mpi_communicator);
-
+
MPI_Op op;
int ierr = MPI_Op_create((MPI_User_function *)&max_reduce, true, &op);
AssertThrow(ierr == MPI_SUCCESS, ExcInternalError());
-
+
MinMaxAvg in;
in.sum = in.min = in.max = my_value;
in.min_index = in.max_index = my_id;
}
#else
-
+
unsigned int n_mpi_processes (const MPI_Comm &)
{
return 1;
}
-
-
-
+
+
+
unsigned int this_mpi_process (const MPI_Comm &)
{
return 0;
}
-
-
+
+
MPI_Comm duplicate_communicator (const MPI_Comm &mpi_communicator)
{
return mpi_communicator;
}
-
-
-
-
+
+
+
+
MinMaxAvg
min_max_avg(const double my_value,
const MPI_Comm &)
{
MinMaxAvg result;
-
+
result.sum = my_value;
result.avg = my_value;
result.min = my_value;
MPI_Initialized(&MPI_has_been_started);
AssertThrow (MPI_has_been_started == 0,
ExcMessage ("MPI error. You can only start MPI once!"));
-
+
int mpi_err;
mpi_err = MPI_Init (&argc, &argv);
AssertThrow (mpi_err == 0,
(void)argv;
(void)owns_mpi;
#endif
-
+
constructor_has_already_run = true;
}
-
-
+
+
MPI_InitFinalize::~MPI_InitFinalize()
{
#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-
+
// make memory pool release all MPI-based vectors that are no
// longer used at this point. this is relevant because the
// static object destructors run for these vectors at the end of
// to errors
# if defined(DEAL_II_USE_TRILINOS)
GrowingVectorMemory<TrilinosWrappers::MPI::Vector>
- ::release_unused_memory ();
+ ::release_unused_memory ();
GrowingVectorMemory<TrilinosWrappers::MPI::BlockVector>
- ::release_unused_memory ();
+ ::release_unused_memory ();
# endif
-
+
// Same for PETSc. only do this if PETSc hasn't been
// terminated yet since PETSc deletes all vectors that
// have been allocated but not freed at the time of
// yield errors of double deallocations
#ifdef DEAL_II_USE_PETSC
if ((PetscInitializeCalled == PETSC_TRUE)
- &&
- (PetscFinalizeCalled == PETSC_FALSE))
- {
- GrowingVectorMemory<PETScWrappers::MPI::Vector>
- ::release_unused_memory ();
- GrowingVectorMemory<PETScWrappers::MPI::BlockVector>
- ::release_unused_memory ();
- GrowingVectorMemory<PETScWrappers::Vector>
- ::release_unused_memory ();
- GrowingVectorMemory<PETScWrappers::BlockVector>
- ::release_unused_memory ();
-
+ &&
+ (PetscFinalizeCalled == PETSC_FALSE))
+ {
+ GrowingVectorMemory<PETScWrappers::MPI::Vector>
+ ::release_unused_memory ();
+ GrowingVectorMemory<PETScWrappers::MPI::BlockVector>
+ ::release_unused_memory ();
+ GrowingVectorMemory<PETScWrappers::Vector>
+ ::release_unused_memory ();
+ GrowingVectorMemory<PETScWrappers::BlockVector>
+ ::release_unused_memory ();
+
# ifdef DEAL_II_USE_SLEPC
- // and now end SLEPc (with PETSc)
- SlepcFinalize();
+ // and now end SLEPc (with PETSc)
+ SlepcFinalize();
# else
- // or just end PETSc.
- PetscFinalize();
+ // or just end PETSc.
+ PetscFinalize();
# endif
- }
+ }
#else
-
+
int mpi_err = 0;
-
+
int MPI_has_been_started = 0;
MPI_Initialized(&MPI_has_been_started);
if (Utilities::System::job_supports_mpi() == true && owns_mpi == true &&
else
mpi_err = MPI_Finalize();
}
-
-
+
+
AssertThrow (mpi_err == 0,
ExcMessage ("An error occurred while calling MPI_Finalize()"));
#endif
#endif
}
-
-
+
+
} // end of namespace MPI
-
+
} // end of namespace Utilities
DEAL_II_NAMESPACE_CLOSE
(digits==3 && i>=1000) ||
(digits==4 && i>=10000)||
(digits==5 && i>=100000)||
- (digits==6 && i>=1000000)||
+ (digits==6 && i>=1000000)||
(digits==7 && i>=10000000)||
(digits==8 && i>=100000000)||
- (digits==9 && i>=1000000000)||
+ (digits==9 && i>=1000000000)||
(i>=1000000000)),
ExcInvalidNumber2StringConversersion(i, digits));
Assert (tr != 0, ExcInternalError());
AssertThrow(
- (tr->settings & parallel::distributed::Triangulation< dim, spacedim >::construct_multigrid_hierarchy),
- ExcMessage("Multigrid DoFs can only be distributed on a parallel Triangulation if the flag construct_multigrid_hierarchy is set in the constructor."));
+ (tr->settings & parallel::distributed::Triangulation< dim, spacedim >::construct_multigrid_hierarchy),
+ ExcMessage("Multigrid DoFs can only be distributed on a parallel Triangulation if the flag construct_multigrid_hierarchy is set in the constructor."));
const unsigned int
// parallel::distributed::Triangulation in serial
{
typedef parallel::distributed::Triangulation<DH::dimension,DH::space_dimension> PTRIA;
- const PTRIA* ptria_p = dynamic_cast<const PTRIA*> (&dof_handler.get_tria());
+ const PTRIA *ptria_p = dynamic_cast<const PTRIA *> (&dof_handler.get_tria());
Assert ((ptria_p == 0 || Utilities::MPI::n_mpi_processes(ptria_p->get_communicator()) == 1),
- ExcMessage ("This function can not be used with distributed triangulations."
- "See the documentation for more information."));
+ ExcMessage ("This function can not be used with distributed triangulations."
+ "See the documentation for more information."));
}
#endif
// parallel::distributed::Triangulation in serial
{
typedef typename parallel::distributed::Triangulation<DH::dimension,DH::space_dimension> PTRIA;
- const PTRIA* ptria_p = dynamic_cast<const PTRIA*> (&dof_handler.get_tria());
+ const PTRIA *ptria_p = dynamic_cast<const PTRIA *> (&dof_handler.get_tria());
Assert ((ptria_p == 0 || Utilities::MPI::n_mpi_processes(ptria_p->get_communicator()) == 1),
- ExcMessage ("This function can not be used with distributed triangulations."
- "See the documentation for more information."));
+ ExcMessage ("This function can not be used with distributed triangulations."
+ "See the documentation for more information."));
}
#endif
for (unsigned int d=0; d<dim; ++d)
indices[d] = 0;
}
-
+
// in initialize_embedding() and
// initialize_restriction(), want to undo
// tensorization on inner loops for
++indices[0];
for (int d=0; d<dim-1; ++d)
if (indices[d]==dofs1d)
- {
- indices[d] = 0;
- indices[d+1]++;
- }
+ {
+ indices[d] = 0;
+ indices[d+1]++;
+ }
}
-
+
// in initialize_embedding() and
// initialize_restriction(), want to undo
// tensorization on inner loops for
unit_support_points.size());
std::vector<Point<1> > points1d (dofs1d);
for (unsigned int i=0; i<dofs1d; ++i)
- {
- const unsigned int j = index_map_inverse[i];
- points1d[i] = Point<1>(unit_support_points[j](0));
- for (unsigned int d=1; d<dim; ++d)
- Assert (unit_support_points[j][d] == 0.,
- ExcInternalError());
- }
+ {
+ const unsigned int j = index_map_inverse[i];
+ points1d[i] = Point<1>(unit_support_points[j](0));
+ for (unsigned int d=1; d<dim; ++d)
+ Assert (unit_support_points[j][d] == 0.,
+ ExcInternalError());
+ }
return Polynomials::generate_complete_Lagrange_basis (points1d);
}
}
FiniteElementData<dim>::L2),
get_riaf_vector(degree),
std::vector<ComponentMask>(1ll , std::vector<bool>(1,true))),
- face_index_map(FE_Q_DG0_Helper::invert_numbering(
- face_lexicographic_to_hierarchic_numbering (degree)))
+ face_index_map(FE_Q_DG0_Helper::invert_numbering(
+ face_lexicographic_to_hierarchic_numbering (degree)))
{
Assert (degree > 0,
ExcMessage ("This element can only be used for polynomial degrees "
this->reinit_restriction_and_prolongation_matrices();
initialize_embedding();
initialize_restriction();
-
+
initialize_quad_dof_index_permutation();
}
get_riaf_vector(points.size()-1),
std::vector<ComponentMask>(1, std::vector<bool>(1,true))),
- face_index_map(FE_Q_DG0_Helper::invert_numbering(
- face_lexicographic_to_hierarchic_numbering (points.size()-1)))
+ face_index_map(FE_Q_DG0_Helper::invert_numbering(
+ face_lexicographic_to_hierarchic_numbering (points.size()-1)))
{
const int degree = points.size()-1;
initialize_embedding();
initialize_restriction();
- initialize_quad_dof_index_permutation();
+ initialize_quad_dof_index_permutation();
}
{
Assert (values.size() == this->unit_support_points.size(),
ExcDimensionMismatch(values.size(),
- this->unit_support_points.size()));
+ this->unit_support_points.size()));
Assert (local_dofs.size() == this->dofs_per_cell,
ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
Assert (this->n_components() == 1,
ExcDimensionMismatch(this->n_components(), 1));
-
+
std::copy(values.begin(), values.end(), local_dofs.begin());
//We don't need the discontinuous function for local interpolation
local_dofs[local_dofs.size()-1]=0.;
{
Assert (values.size() == this->unit_support_points.size(),
ExcDimensionMismatch(values.size(),
- this->unit_support_points.size()));
+ this->unit_support_points.size()));
Assert (local_dofs.size() == this->dofs_per_cell,
ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
Assert (values[0].size() >= offset+this->n_components(),
ExcDimensionMismatch(values[0].size(),offset+this->n_components()));
-
+
for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
- {
- const std::pair<unsigned int, unsigned int> index
- = this->system_to_component_index(i);
- local_dofs[i] = values[i](offset+index.first);
- }
+ {
+ const std::pair<unsigned int, unsigned int> index
+ = this->system_to_component_index(i);
+ local_dofs[i] = values[i](offset+index.first);
+ }
//We don't need the discontinuous function for local interpolation
local_dofs[local_dofs.size()-1]=0.;
}
-
-
-
-
+
+
+
+
template <int dim, int spacedim>
void
FE_Q_DG0<dim,spacedim>::interpolate(
{
Assert (values[0].size() == this->unit_support_points.size(),
ExcDimensionMismatch(values.size(),
- this->unit_support_points.size()));
+ this->unit_support_points.size()));
Assert (local_dofs.size() == this->dofs_per_cell,
ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
Assert (values.size() == this->n_components(),
ExcDimensionMismatch(values.size(), this->n_components()));
-
+
for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
- {
- const std::pair<unsigned int, unsigned int> index
- = this->system_to_component_index(i);
- local_dofs[i] = values[index.first][i];
- }
+ {
+ const std::pair<unsigned int, unsigned int> index
+ = this->system_to_component_index(i);
+ local_dofs[i] = values[index.first][i];
+ }
//We don't need the discontinuous function for local interpolation
local_dofs[local_dofs.size()-1]=0.;
}
const Point<dim> p = this->unit_support_points[j];
for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
cell_interpolation(j,i) = this->poly_space.compute_value (i, p);
-
+
for (unsigned int i=0; i<source_fe.dofs_per_cell-1; ++i)
source_interpolation(j,i) = source_fe.poly_space.compute_value (i,p);
-
+
}
- //the discontinuous node is not transformed
- cell_interpolation(this->dofs_per_cell-1,
- this->dofs_per_cell-1)=1.;
- source_interpolation(this->dofs_per_cell-1,
- source_fe.dofs_per_cell-1)=1.;
+ //the discontinuous node is not transformed
+ cell_interpolation(this->dofs_per_cell-1,
+ this->dofs_per_cell-1)=1.;
+ source_interpolation(this->dofs_per_cell-1,
+ source_fe.dofs_per_cell-1)=1.;
// then compute the
// interpolation matrix
// FE_Q_DG0s or if the other one is an
// FE_Nothing
if (const FE_Q_DG0<dim,spacedim> *fe_q_dg0_other = dynamic_cast<const
-FE_Q_DG0<dim,spacedim>*>(&fe_other))
+ FE_Q_DG0<dim,spacedim> *>(&fe_other))
{
// dofs are located along lines, so two
// dofs are identical if they are
for (unsigned int j=0; j<q-1; ++j)
if (std::fabs(this->unit_support_points[index_map_inverse[i+1]][0]-
-fe_q_dg0_other->unit_support_points[index_map_inverse_other[j+1]][0])
+ fe_q_dg0_other->unit_support_points[index_map_inverse_other[j+1]][0])
< 1e-14)
identities.push_back (std::make_pair(i,j));
// FE_Q_DG0s or if the other one is an
// FE_Nothing
if (const FE_Q_DG0<dim,spacedim> *fe_q_dg0_other = dynamic_cast<const
-FE_Q_DG0<dim,spacedim>*>(&fe_other))
+ FE_Q_DG0<dim,spacedim> *>(&fe_other))
{
// this works exactly like the line
// case above, except that now we have
for (unsigned int j2=0; j2<q-1; ++j2)
if ((std::fabs(this->unit_support_points[index_map_inverse[i1+1]][0]-
-fe_q_dg0_other->unit_support_points[index_map_inverse_other[j1+1]][0])
+ fe_q_dg0_other->unit_support_points[index_map_inverse_other[j1+1]][0])
< 1e-14)
&&
(std::fabs(this->unit_support_points[index_map_inverse[i2+1]][0]-
-fe_q_dg0_other->unit_support_points[index_map_inverse_other[j2+1]][0])
+ fe_q_dg0_other->unit_support_points[index_map_inverse_other[j2+1]][0])
< 1e-14))
identities.push_back (std::make_pair(i1*(p-1)+i2,
j1*(q-1)+j2));
FE_Q_DG0<dim,spacedim>::get_riaf_vector(const unsigned int deg)
{
std::vector<bool> riaf
- (FiniteElementData<dim> (get_dpo_vector(deg),1,deg).dofs_per_cell,
- false);
+ (FiniteElementData<dim> (get_dpo_vector(deg),1,deg).dofs_per_cell,
+ false);
riaf[riaf.size()-1]=true;
return riaf;
}
// cell matrix. The value eps is used a threshold to decide when certain
// evaluations of the Lagrange polynomials are zero or one.
const double eps = 1e-15*this->degree*dim;
-
- #ifdef DEBUG
+
+#ifdef DEBUG
// in DEBUG mode, check that the evaluation of support points (except for
// the discontinuous node) in the current numbering gives the identity
// operation
for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
- {
- Assert (std::fabs (1.-this->poly_space.compute_value
- (i, this->unit_support_points[i])) < eps,
- ExcInternalError());
- for (unsigned int j=0; j<this->dofs_per_cell-1; ++j)
- if (j!=i)
- Assert (std::fabs (this->poly_space.compute_value
- (i, this->unit_support_points[j])) < eps,
- ExcInternalError());
- }
- #endif
-
+ {
+ Assert (std::fabs (1.-this->poly_space.compute_value
+ (i, this->unit_support_points[i])) < eps,
+ ExcInternalError());
+ for (unsigned int j=0; j<this->dofs_per_cell-1; ++j)
+ if (j!=i)
+ Assert (std::fabs (this->poly_space.compute_value
+ (i, this->unit_support_points[j])) < eps,
+ ExcInternalError());
+ }
+#endif
+
// to efficiently evaluate the polynomial at the subcell, make use of the
// tensor product structure of this element and only evaluate 1D
// information from the polynomial. This makes the cost of this function
// almost negligible also for high order elements
const unsigned int dofs1d = this->degree+1;
std::vector<Table<2,double> > subcell_evaluations
- (dim, Table<2,double>(dofs1d, dofs1d));
+ (dim, Table<2,double>(dofs1d, dofs1d));
const std::vector<unsigned int> &index_map_inverse =
this->poly_space.get_numbering_inverse();
-
+
// recreate 1D polynomials, the last entry is not used
std::vector<Polynomials::Polynomial<double> > poly_space1d =
FE_Q_DG0_Helper::generate_poly_space1d (this->unit_support_points,
- index_map_inverse, dofs1d);
-
+ index_map_inverse, dofs1d);
+
// helper value: step size how to walk through diagonal and how many
// points we have left apart from the first dimension
unsigned int step_size_diag = 0;
{
unsigned int factor = 1;
for (unsigned int d=0; d<dim; ++d)
- {
- step_size_diag += factor;
- factor *= dofs1d;
- }
+ {
+ step_size_diag += factor;
+ factor *= dofs1d;
+ }
}
-
+
// next evaluate the functions for the different refinement cases.
for (unsigned int ref=0; ref<RefinementCase<dim>::isotropic_refinement; ++ref)
for (unsigned int child=0;
child<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref+1));
++child)
- {
- // go through the points in diagonal to capture variation in all
- // directions simultaneously
- for (unsigned int j=0; j<dofs1d; ++j)
- {
- const unsigned int diag_comp = index_map_inverse[j*step_size_diag];
- const Point<dim> p_subcell = this->unit_support_points[diag_comp];
- const Point<dim> p_cell =
+ {
+ // go through the points in diagonal to capture variation in all
+ // directions simultaneously
+ for (unsigned int j=0; j<dofs1d; ++j)
+ {
+ const unsigned int diag_comp = index_map_inverse[j*step_size_diag];
+ const Point<dim> p_subcell = this->unit_support_points[diag_comp];
+ const Point<dim> p_cell =
GeometryInfo<dim>::child_to_cell_coordinates
- (p_subcell, child, RefinementCase<dim>(ref+1));
- for (unsigned int i=0; i<dofs1d; ++i)
- for (unsigned int d=0; d<dim; ++d)
- {
- const double cell_value = poly_space1d[i].value (p_cell[d]);
-
- // cut off values that are too small. note that we have here
- // Lagrange interpolation functions, so they should be zero at
- // almost all points, and one at the others, at least on the
- // subcells. so set them to their exact values.
- //
- // the actual cut-off value is somewhat fuzzy, but it works for
- // 2e-13*degree*dim (see above), which is kind of reasonable
- // given that we compute the values of the polynomials via an
- // degree-step recursion and then multiply the 1d-values. this
- // gives us a linear growth in degree*dim, times a small
- // constant.
- //
- // the embedding matrix is given by applying the inverse of the
- // subcell matrix on the cell_interpolation matrix. since the
- // subcell matrix is actually only a permutation vector, all we
- // need to do is to switch the rows we write the data into.
- // moreover, cut off very small values here
- if (std::fabs(cell_value) < eps)
- subcell_evaluations[d](j,i) = 0;
- else
- subcell_evaluations[d](j,i) = cell_value;
- }
- }
-
- // now expand from 1D info. block innermost dimension (x_0) in order
- // to avoid difficult checks at innermost loop
- unsigned int j_indices[dim];
- FE_Q_DG0_Helper::zero_indices<dim> (j_indices);
- for (unsigned int j=0; j<this->dofs_per_cell-1; j+=dofs1d)
- {
- unsigned int i_indices[dim];
- FE_Q_DG0_Helper::zero_indices<dim> (i_indices);
- for (unsigned int i=0; i<this->dofs_per_cell-1; i+=dofs1d)
- {
- double val_extra_dim = 1.;
- for (unsigned int d=1; d<dim; ++d)
- val_extra_dim *= subcell_evaluations[d](j_indices[d-1],
- i_indices[d-1]);
-
- // innermost sum where we actually compute. the same as
- // this->prolongation[ref][child](j,i) =
- // this->poly_space.compute_value (i, p_cell);
- for (unsigned int jj=0; jj<dofs1d; ++jj)
- {
- const unsigned int j_ind = index_map_inverse[j+jj];
- for (unsigned int ii=0; ii<dofs1d; ++ii)
-
- this->prolongation[ref][child](j_ind,index_map_inverse[i+ii])
- = val_extra_dim * subcell_evaluations[0](jj,ii);
- }
-
- // update indices that denote the tensor product position. a bit
- // fuzzy and therefore not done for innermost x_0 direction
- FE_Q_DG0_Helper::increment_indices<dim> (i_indices, dofs1d);
- }
- Assert (i_indices[dim-1] == 1, ExcInternalError());
- FE_Q_DG0_Helper::increment_indices<dim> (j_indices, dofs1d);
- }
-
- // the discontinuous node is simply mapped on the discontinuous node
- // on the child cell
- for (unsigned int j=0; j<this->dofs_per_cell-1; j++)
- {
- this->prolongation[ref][child](j,this->dofs_per_cell-1) = 0.;
- this->prolongation[ref][child](this->dofs_per_cell-1,j) = 0.;
- }
- this->prolongation[ref][child](this->dofs_per_cell-1,
- this->dofs_per_cell-1) = 1.;
-
- // make sure that the row sum is one. this must be so since for this
- // element, the continuous shape functions add up to one and the
- // discontinuous node is mapped to the discontinuous node on the
- // child cell
- #ifdef DEBUG
- for (unsigned int row=0; row<this->dofs_per_cell; ++row)
- {
- double sum = 0;
- for (unsigned int col=0; col<this->dofs_per_cell; ++col)
- sum += this->prolongation[ref][child](row,col);
- Assert (std::fabs(sum-1.) < eps, ExcInternalError());
- }
- #endif
- }
+ (p_subcell, child, RefinementCase<dim>(ref+1));
+ for (unsigned int i=0; i<dofs1d; ++i)
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ const double cell_value = poly_space1d[i].value (p_cell[d]);
+
+ // cut off values that are too small. note that we have here
+ // Lagrange interpolation functions, so they should be zero at
+ // almost all points, and one at the others, at least on the
+ // subcells. so set them to their exact values.
+ //
+ // the actual cut-off value is somewhat fuzzy, but it works for
+ // 2e-13*degree*dim (see above), which is kind of reasonable
+ // given that we compute the values of the polynomials via an
+ // degree-step recursion and then multiply the 1d-values. this
+ // gives us a linear growth in degree*dim, times a small
+ // constant.
+ //
+ // the embedding matrix is given by applying the inverse of the
+ // subcell matrix on the cell_interpolation matrix. since the
+ // subcell matrix is actually only a permutation vector, all we
+ // need to do is to switch the rows we write the data into.
+ // moreover, cut off very small values here
+ if (std::fabs(cell_value) < eps)
+ subcell_evaluations[d](j,i) = 0;
+ else
+ subcell_evaluations[d](j,i) = cell_value;
+ }
+ }
+
+ // now expand from 1D info. block innermost dimension (x_0) in order
+ // to avoid difficult checks at innermost loop
+ unsigned int j_indices[dim];
+ FE_Q_DG0_Helper::zero_indices<dim> (j_indices);
+ for (unsigned int j=0; j<this->dofs_per_cell-1; j+=dofs1d)
+ {
+ unsigned int i_indices[dim];
+ FE_Q_DG0_Helper::zero_indices<dim> (i_indices);
+ for (unsigned int i=0; i<this->dofs_per_cell-1; i+=dofs1d)
+ {
+ double val_extra_dim = 1.;
+ for (unsigned int d=1; d<dim; ++d)
+ val_extra_dim *= subcell_evaluations[d](j_indices[d-1],
+ i_indices[d-1]);
+
+ // innermost sum where we actually compute. the same as
+ // this->prolongation[ref][child](j,i) =
+ // this->poly_space.compute_value (i, p_cell);
+ for (unsigned int jj=0; jj<dofs1d; ++jj)
+ {
+ const unsigned int j_ind = index_map_inverse[j+jj];
+ for (unsigned int ii=0; ii<dofs1d; ++ii)
+
+ this->prolongation[ref][child](j_ind,index_map_inverse[i+ii])
+ = val_extra_dim * subcell_evaluations[0](jj,ii);
+ }
+
+ // update indices that denote the tensor product position. a bit
+ // fuzzy and therefore not done for innermost x_0 direction
+ FE_Q_DG0_Helper::increment_indices<dim> (i_indices, dofs1d);
+ }
+ Assert (i_indices[dim-1] == 1, ExcInternalError());
+ FE_Q_DG0_Helper::increment_indices<dim> (j_indices, dofs1d);
+ }
+
+ // the discontinuous node is simply mapped on the discontinuous node
+ // on the child cell
+ for (unsigned int j=0; j<this->dofs_per_cell-1; j++)
+ {
+ this->prolongation[ref][child](j,this->dofs_per_cell-1) = 0.;
+ this->prolongation[ref][child](this->dofs_per_cell-1,j) = 0.;
+ }
+ this->prolongation[ref][child](this->dofs_per_cell-1,
+ this->dofs_per_cell-1) = 1.;
+
+ // make sure that the row sum is one. this must be so since for this
+ // element, the continuous shape functions add up to one and the
+ // discontinuous node is mapped to the discontinuous node on the
+ // child cell
+#ifdef DEBUG
+ for (unsigned int row=0; row<this->dofs_per_cell; ++row)
+ {
+ double sum = 0;
+ for (unsigned int col=0; col<this->dofs_per_cell; ++col)
+ sum += this->prolongation[ref][child](row,col);
+ Assert (std::fabs(sum-1.) < eps, ExcInternalError());
+ }
+#endif
+ }
}
// about this.
// for the discontinuous node we just take the mean of all the child cells'
// contributions.
-
+
const double eps = 1e-15*this->degree*dim;
const std::vector<unsigned int> &index_map_inverse =
- this->poly_space.get_numbering_inverse();
-
+ this->poly_space.get_numbering_inverse();
+
// recreate 1D polynomials for faster evaluation of polynomial, the last
// entry is not used
const unsigned int dofs1d = this->degree+1;
FE_Q_DG0_Helper::generate_poly_space1d (this->unit_support_points,
index_map_inverse, dofs1d);
std::vector<Tensor<1,dim> > evaluations1d (dofs1d);
-
+
for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
- {
- unsigned int mother_dof = index_map_inverse[i];
- const Point<dim> p_cell = this->unit_support_points[mother_dof];
-
- // then find the children on which the interpolation point is located
- for (unsigned int ref=RefinementCase<dim>::cut_x;
- ref<=RefinementCase<dim>::isotropic_refinement; ++ref)
- for (unsigned int child=0;
- child<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref));
- ++child)
+ {
+ unsigned int mother_dof = index_map_inverse[i];
+ const Point<dim> p_cell = this->unit_support_points[mother_dof];
+
+ // then find the children on which the interpolation point is located
+ for (unsigned int ref=RefinementCase<dim>::cut_x;
+ ref<=RefinementCase<dim>::isotropic_refinement; ++ref)
+ for (unsigned int child=0;
+ child<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref));
+ ++child)
+ {
+ // check whether this interpolation point is inside this child
+ // cell
+ const Point<dim> p_subcell =
+ GeometryInfo<dim>::cell_to_child_coordinates
+ (p_cell, child, RefinementCase<dim>(ref));
+ if (GeometryInfo<dim>::is_inside_unit_cell (p_subcell))
{
- // check whether this interpolation point is inside this child
- // cell
- const Point<dim> p_subcell =
- GeometryInfo<dim>::cell_to_child_coordinates
- (p_cell, child, RefinementCase<dim>(ref));
- if (GeometryInfo<dim>::is_inside_unit_cell (p_subcell))
- {
- // same logic as in initialize_embedding to evaluate the
- // polynomial faster than from the tensor product: since we
- // evaluate all polynomials, it is much faster to just
- // compute the 1D values for all polynomials before and then
- // get the dim-data.
- for (unsigned int j=0; j<dofs1d; ++j)
- for (unsigned int d=0; d<dim; ++d)
- evaluations1d[j][d] = poly_space1d[j].value(p_subcell[d]);
- unsigned int j_indices[dim];
- FE_Q_DG0_Helper::zero_indices<dim> (j_indices);
- double sum_check = 0;
- for (unsigned int j = 0; j<this->dofs_per_cell-1; j += dofs1d)
+ // same logic as in initialize_embedding to evaluate the
+ // polynomial faster than from the tensor product: since we
+ // evaluate all polynomials, it is much faster to just
+ // compute the 1D values for all polynomials before and then
+ // get the dim-data.
+ for (unsigned int j=0; j<dofs1d; ++j)
+ for (unsigned int d=0; d<dim; ++d)
+ evaluations1d[j][d] = poly_space1d[j].value(p_subcell[d]);
+ unsigned int j_indices[dim];
+ FE_Q_DG0_Helper::zero_indices<dim> (j_indices);
+ double sum_check = 0;
+ for (unsigned int j = 0; j<this->dofs_per_cell-1; j += dofs1d)
{
double val_extra_dim = 1.;
for (unsigned int d=1; d<dim; ++d)
val_extra_dim *= evaluations1d[j_indices[d-1]][d];
for (unsigned int jj=0; jj<dofs1d; ++jj)
- {
-
- // find the child shape function(s) corresponding
- // to this point. Usually this is just one function;
- // however, when we use FE_Q_DG0 on arbitrary nodes a
- // parent support point might not be a child support
- // point, and then we will get more than one nonzero
- // value per row. Still, the values should sum up to 1
- const double val = val_extra_dim * evaluations1d[jj][0];
- const unsigned int child_dof = index_map_inverse[j+jj];
- if (std::fabs (val-1.) < eps)
- this->restriction[ref-1][child]
- (mother_dof,child_dof)=1.;
- else if (std::fabs(val) > eps)
- this->restriction[ref-1][child]
- (mother_dof,child_dof)=val;
- sum_check += val;
- }
+ {
+
+ // find the child shape function(s) corresponding
+ // to this point. Usually this is just one function;
+ // however, when we use FE_Q_DG0 on arbitrary nodes a
+ // parent support point might not be a child support
+ // point, and then we will get more than one nonzero
+ // value per row. Still, the values should sum up to 1
+ const double val = val_extra_dim * evaluations1d[jj][0];
+ const unsigned int child_dof = index_map_inverse[j+jj];
+ if (std::fabs (val-1.) < eps)
+ this->restriction[ref-1][child]
+ (mother_dof,child_dof)=1.;
+ else if (std::fabs(val) > eps)
+ this->restriction[ref-1][child]
+ (mother_dof,child_dof)=val;
+ sum_check += val;
+ }
FE_Q_DG0_Helper::increment_indices<dim> (j_indices, dofs1d);
}
- Assert (std::fabs(sum_check-1) < eps,
- ExcInternalError());
- }
+ Assert (std::fabs(sum_check-1) < eps,
+ ExcInternalError());
}
- }
+ }
+ }
// for the discontinuous node we take the mean of all the child cells'
// contributions. so there is also a diagonal entry equal to the inverse of
// n_children.
for (unsigned int ref=RefinementCase<dim>::cut_x;
ref<=RefinementCase<dim>::isotropic_refinement; ++ref)
- for (unsigned int child=0;
- child<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref));
- ++child)
- {
- for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
- {
- this->restriction[ref-1][child](i,this->dofs_per_cell-1)=0.;
- this->restriction[ref-1][child](this->dofs_per_cell-1,i)=0.;
- }
- this->restriction[ref-1][child]
- (this->dofs_per_cell-1,this->dofs_per_cell-1) =
- 1./(double)GeometryInfo<dim>::n_children
- (RefinementCase<dim>(ref));
- }
-
+ for (unsigned int child=0;
+ child<GeometryInfo<dim>::n_children(RefinementCase<dim>(ref));
+ ++child)
+ {
+ for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
+ {
+ this->restriction[ref-1][child](i,this->dofs_per_cell-1)=0.;
+ this->restriction[ref-1][child](this->dofs_per_cell-1,i)=0.;
+ }
+ this->restriction[ref-1][child]
+ (this->dofs_per_cell-1,this->dofs_per_cell-1) =
+ 1./(double)GeometryInfo<dim>::n_children
+ (RefinementCase<dim>(ref));
+ }
+
}
-
+
//---------------------------------------------------------------------------
template<>
void
GridGenerator::parallelogram (Triangulation<2> &tria,
- const Point<2> (&corners)[2],
- const bool colorize)
+ const Point<2> (&corners)[2],
+ const bool colorize)
{
std::vector<Point<2> > vertices (GeometryInfo<2>::vertices_per_cell);
template<>
void
GridGenerator::parallelogram (Triangulation<2> &tria,
- const Tensor<2,2> &corners,
- const bool colorize)
+ const Tensor<2,2> &corners,
+ const bool colorize)
{
// simply pass everything to the other function of same name
const Point<2> x[2] = { corners[0], corners[1] };
template<int dim>
void
GridGenerator::parallelepiped (Triangulation<dim> &tria,
- const Point<dim> (&corners) [dim],
- const bool colorize)
+ const Point<dim> (&corners) [dim],
+ const bool colorize)
{
// Check that none of the user defined vertices overlap
for (unsigned int i=0; i<dim; ++i)
template<int dim>
void
GridGenerator::subdivided_parallelepiped (Triangulation<dim> &tria,
- const unsigned int n_subdivisions,
- const Point<dim> (&corners) [dim],
- const bool colorize)
+ const unsigned int n_subdivisions,
+ const Point<dim> (&corners) [dim],
+ const bool colorize)
{
// Equalise number of subdivisions in each dim-direction, heir
// validity will be checked later
n_subdivisions_[i] = n_subdivisions;
// and call the function below
- GridGenerator::subdivided_parallelepiped (tria, n_subdivisions_,
- corners,
- colorize);
+ GridGenerator::subdivided_parallelepiped (tria, n_subdivisions_,
+ corners,
+ colorize);
}
template<int dim>
void
GridGenerator::subdivided_parallelepiped (Triangulation<dim> &tria,
- const unsigned int ( n_subdivisions) [dim],
- const Point<dim> (&corners) [dim],
- const bool colorize)
+ const unsigned int ( n_subdivisions) [dim],
+ const Point<dim> (&corners) [dim],
+ const bool colorize)
{
// Zero n_subdivisions is the origin only, which makes no sense
for (unsigned int i=0; i<dim; ++i)
{
case 1:
for (unsigned int x=0; x<=n_subdivisions[0]; ++x)
- points.push_back (Point<dim> (x*delta[0]));
+ points.push_back (Point<dim> (x*delta[0]));
break;
-
+
case 2:
for (unsigned int y=0; y<=n_subdivisions[1]; ++y)
- for (unsigned int x=0; x<=n_subdivisions[0]; ++x)
- points.push_back (Point<dim> (x*delta[0] + y*delta[1]));
+ for (unsigned int x=0; x<=n_subdivisions[0]; ++x)
+ points.push_back (Point<dim> (x*delta[0] + y*delta[1]));
break;
case 3:
for (unsigned int z=0; z<=n_subdivisions[2]; ++z)
- for (unsigned int y=0; y<=n_subdivisions[1]; ++y)
- for (unsigned int x=0; x<=n_subdivisions[0]; ++x)
- points.push_back (Point<dim> (x*delta[0] + y*delta[1] + z*delta[2]));
+ for (unsigned int y=0; y<=n_subdivisions[1]; ++y)
+ for (unsigned int x=0; x<=n_subdivisions[0]; ++x)
+ points.push_back (Point<dim> (x*delta[0] + y*delta[1] + z*delta[2]));
break;
default:
for (unsigned int i=0; i<dim; ++i)
n_points *= n_subdivisions[i]+1;
- points.erase (std::unique (points.begin (), points.end ()),
- points.end ());
- Assert (points.size ()==n_points,
- ExcInternalError());
+ points.erase (std::unique (points.begin (), points.end ()),
+ points.end ());
+ Assert (points.size ()==n_points,
+ ExcInternalError());
#endif
-
- // Prepare cell data
+
+ // Prepare cell data
unsigned int n_cells = 1;
for (unsigned int i=0; i<dim; ++i)
n_cells *= n_subdivisions[i];
std::vector<CellData<dim> > cells (n_cells);
- // Create fixed ordering of
+ // Create fixed ordering of
switch (dim)
{
case 1:
cells[x].vertices[0] = x;
cells[x].vertices[1] = x+1;
- // wipe material id
+ // wipe material id
cells[x].material_id = 0;
}
break;
case 2:
- {
- // Shorthand
- const unsigned int n_dy = n_subdivisions[1];
- const unsigned int n_dx = n_subdivisions[0];
-
- for (unsigned int y=0; y<n_dy; ++y)
- for (unsigned int x=0; x<n_dx; ++x)
- {
- const unsigned int c = y*n_dx + x;
- cells[c].vertices[0] = y*(n_dx+1) + x;
- cells[c].vertices[1] = y*(n_dx+1) + x+1;
- cells[c].vertices[2] = (y+1)*(n_dx+1) + x;
- cells[c].vertices[3] = (y+1)*(n_dx+1) + x+1;
-
- // wipe material id
- cells[c].material_id = 0;
- }
- }
- break;
-
+ {
+ // Shorthand
+ const unsigned int n_dy = n_subdivisions[1];
+ const unsigned int n_dx = n_subdivisions[0];
+
+ for (unsigned int y=0; y<n_dy; ++y)
+ for (unsigned int x=0; x<n_dx; ++x)
+ {
+ const unsigned int c = y*n_dx + x;
+ cells[c].vertices[0] = y*(n_dx+1) + x;
+ cells[c].vertices[1] = y*(n_dx+1) + x+1;
+ cells[c].vertices[2] = (y+1)*(n_dx+1) + x;
+ cells[c].vertices[3] = (y+1)*(n_dx+1) + x+1;
+
+ // wipe material id
+ cells[c].material_id = 0;
+ }
+ }
+ break;
+
case 3:
- {
- // Shorthand
- const unsigned int n_dz = n_subdivisions[2];
- const unsigned int n_dy = n_subdivisions[1];
- const unsigned int n_dx = n_subdivisions[0];
-
- for (unsigned int z=0; z<n_dz; ++z)
- for (unsigned int y=0; y<n_dy; ++y)
- for (unsigned int x=0; x<n_dx; ++x)
- {
- const unsigned int c = z*n_dy*n_dx + y*n_dx + x;
-
- cells[c].vertices[0] = z*(n_dy+1)*(n_dx+1) + y*(n_dx+1) + x;
- cells[c].vertices[1] = z*(n_dy+1)*(n_dx+1) + y*(n_dx+1) + x+1;
- cells[c].vertices[2] = z*(n_dy+1)*(n_dx+1) + (y+1)*(n_dx+1) + x;
- cells[c].vertices[3] = z*(n_dy+1)*(n_dx+1) + (y+1)*(n_dx+1) + x+1;
- cells[c].vertices[4] = (z+1)*(n_dy+1)*(n_dx+1) + y*(n_dx+1) + x;
- cells[c].vertices[5] = (z+1)*(n_dy+1)*(n_dx+1) + y*(n_dx+1) + x+1;
- cells[c].vertices[6] = (z+1)*(n_dy+1)*(n_dx+1) + (y+1)*(n_dx+1) + x;
- cells[c].vertices[7] = (z+1)*(n_dy+1)*(n_dx+1) + (y+1)*(n_dx+1) + x+1;
-
- // wipe material id
- cells[c].material_id = 0;
- }
- break;
- }
-
+ {
+ // Shorthand
+ const unsigned int n_dz = n_subdivisions[2];
+ const unsigned int n_dy = n_subdivisions[1];
+ const unsigned int n_dx = n_subdivisions[0];
+
+ for (unsigned int z=0; z<n_dz; ++z)
+ for (unsigned int y=0; y<n_dy; ++y)
+ for (unsigned int x=0; x<n_dx; ++x)
+ {
+ const unsigned int c = z*n_dy*n_dx + y*n_dx + x;
+
+ cells[c].vertices[0] = z*(n_dy+1)*(n_dx+1) + y*(n_dx+1) + x;
+ cells[c].vertices[1] = z*(n_dy+1)*(n_dx+1) + y*(n_dx+1) + x+1;
+ cells[c].vertices[2] = z*(n_dy+1)*(n_dx+1) + (y+1)*(n_dx+1) + x;
+ cells[c].vertices[3] = z*(n_dy+1)*(n_dx+1) + (y+1)*(n_dx+1) + x+1;
+ cells[c].vertices[4] = (z+1)*(n_dy+1)*(n_dx+1) + y*(n_dx+1) + x;
+ cells[c].vertices[5] = (z+1)*(n_dy+1)*(n_dx+1) + y*(n_dx+1) + x+1;
+ cells[c].vertices[6] = (z+1)*(n_dy+1)*(n_dx+1) + (y+1)*(n_dx+1) + x;
+ cells[c].vertices[7] = (z+1)*(n_dy+1)*(n_dx+1) + (y+1)*(n_dx+1) + x+1;
+
+ // wipe material id
+ cells[c].material_id = 0;
+ }
+ break;
+ }
+
default:
Assert (false, ExcNotImplemented());
}
dof_handler.distribute_dofs(q1);
CompressedSparsityPattern c_sparsity_pattern (dof_handler.n_dofs (),
- dof_handler.n_dofs ());
+ dof_handler.n_dofs ());
DoFTools::make_sparsity_pattern (dof_handler, c_sparsity_pattern);
c_sparsity_pattern.compress ();
// Set options
PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
- format);
+ format);
// Write to screen
MatView (matrix, PETSC_VIEWER_STDOUT_WORLD);
out.setf (std::ios::fixed, std::ios::floatfield);
for ( unsigned int i = 0;
- i < Utilities::MPI::n_mpi_processes(communicator);
- i++)
- {
- // This is slow, but most likely only used to debug.
- MPI_Barrier(communicator);
- if (i == Utilities::MPI::this_mpi_process(communicator))
- {
- if (across)
- {
- out << "[Proc" << i << " " << istart << "-" << iend-1 << "]" << ' ';
- for (PetscInt i=0; i<nlocal; ++i)
- out << val[i] << ' ';
- }
- else
- {
- out << "[Proc " << i << " " << istart << "-" << iend-1 << "]" << std::endl;
- for (PetscInt i=0; i<nlocal; ++i)
- out << val[i] << std::endl;
- }
- out << std::endl;
- }
- }
+ i < Utilities::MPI::n_mpi_processes(communicator);
+ i++)
+ {
+ // This is slow, but most likely only used to debug.
+ MPI_Barrier(communicator);
+ if (i == Utilities::MPI::this_mpi_process(communicator))
+ {
+ if (across)
+ {
+ out << "[Proc" << i << " " << istart << "-" << iend-1 << "]" << ' ';
+ for (PetscInt i=0; i<nlocal; ++i)
+ out << val[i] << ' ';
+ }
+ else
+ {
+ out << "[Proc " << i << " " << istart << "-" << iend-1 << "]" << std::endl;
+ for (PetscInt i=0; i<nlocal; ++i)
+ out << val[i] << std::endl;
+ }
+ out << std::endl;
+ }
+ }
// reset output format
out.flags (old_flags);
out.precision(old_precision);
{
int ierr;
- /*
- TODO: PETSc dublicates communicators, so this does not work (you put MPI_COMM_SELF in, but get something other out when you ask PETSc for the communicator. This mainly fails due to the MatrixFree classes, that can not ask PETSc for a communicator. //Timo Heister
+ /*
+ TODO: PETSc dublicates communicators, so this does not work (you put MPI_COMM_SELF in, but get something other out when you ask PETSc for the communicator. This mainly fails due to the MatrixFree classes, that can not ask PETSc for a communicator. //Timo Heister
Assert(A.get_mpi_communicator()==mpi_communicator, ExcMessage("PETSc Solver and Matrix need to use the same MPI_Comm."));
Assert(x.get_mpi_communicator()==mpi_communicator, ExcMessage("PETSc Solver and Vector need to use the same MPI_Comm."));
Assert(b.get_mpi_communicator()==mpi_communicator, ExcMessage("PETSc Solver and Vector need to use the same MPI_Comm."));
- */
+ */
// first create a solver object if this
// is necessary
// Set options
PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
- format);
+ format);
// Write to screen
VecView (vector, PETSC_VIEWER_STDOUT_WORLD);
solver_control (cn),
mpi_communicator (mpi_communicator),
set_which (EPS_LARGEST_MAGNITUDE),
- opA (NULL),
+ opA (NULL),
opB (NULL),
initial_vector (NULL),
transformation (NULL)
// create a solver object if this is necessary
if (solver_data.get() == 0)
{
- // reset solver dtaa
- solver_data.reset (new SolverData());
-
- // create eigensolver context and set operators
- ierr = EPSCreate (mpi_communicator, &solver_data->eps);
- AssertThrow (ierr == 0, ExcSLEPcError(ierr));
-
- // set eigenspectrum problem type (general/standard)
- AssertThrow (opA, ExcSLEPcWrappersUsageError());
- if (opB)
- ierr = EPSSetOperators (solver_data->eps, *opA, *opB);
- else
- ierr = EPSSetOperators (solver_data->eps, *opA, PETSC_NULL);
- AssertThrow (ierr == 0, ExcSLEPcError(ierr));
-
- // set runtime options
- set_solver_type (solver_data->eps);
+ // reset solver dtaa
+ solver_data.reset (new SolverData());
+
+ // create eigensolver context and set operators
+ ierr = EPSCreate (mpi_communicator, &solver_data->eps);
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
+ // set eigenspectrum problem type (general/standard)
+ AssertThrow (opA, ExcSLEPcWrappersUsageError());
+ if (opB)
+ ierr = EPSSetOperators (solver_data->eps, *opA, *opB);
+ else
+ ierr = EPSSetOperators (solver_data->eps, *opA, PETSC_NULL);
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
+ // set runtime options
+ set_solver_type (solver_data->eps);
}
// set the initial vector(s) if any
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
// Set convergence test to be absolute
- ierr = EPSSetConvergenceTest (solver_data->eps, EPS_CONV_ABS);
- AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+ ierr = EPSSetConvergenceTest (solver_data->eps, EPS_CONV_ABS);
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
// Set the convergence test function
- // ierr = EPSSetConvergenceTestFunction (solver_data->eps, &convergence_test,
- // reinterpret_cast<void *>(&solver_control));
- // AssertThrow (ierr == 0, ExcSLEPcError(ierr));
-
+ // ierr = EPSSetConvergenceTestFunction (solver_data->eps, &convergence_test,
+ // reinterpret_cast<void *>(&solver_control));
+ // AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
// solve the eigensystem
ierr = EPSSolve (solver_data->eps);
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
#ifdef PETSC_USE_64BIT_INDICES
reinterpret_cast<PetscInt *>(n_converged)
#else
- reinterpret_cast<int *>(n_converged)
+ reinterpret_cast<int *>(n_converged)
#endif
- );
+ );
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
int n_iterations = 0;
// complete eigenspectrum
{
// get the number of solver iterations
- ierr = EPSGetIterationNumber (solver_data->eps, &n_iterations);
- AssertThrow (ierr == 0, ExcSLEPcError(ierr));
-
+ ierr = EPSGetIterationNumber (solver_data->eps, &n_iterations);
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+
// get the residual norm of the most extreme eigenvalue
ierr = EPSComputeResidualNorm (solver_data->eps, 0, &residual_norm);
- AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
// check the solver state
const SolverControl::State state
- = solver_control.check (n_iterations, residual_norm);
+ = solver_control.check (n_iterations, residual_norm);
// get the solver state according to SLEPc
get_solver_state (state);
-
+
// and in case of failure: throw exception
- if (solver_control.last_check () != SolverControl::success)
- throw SolverControl::NoConvergence (solver_control.last_step (),
- solver_control.last_value ());
+ if (solver_control.last_check () != SolverControl::success)
+ throw SolverControl::NoConvergence (solver_control.last_step (),
+ solver_control.last_value ());
}
}
{
if (solver_data.get () == 0)
return NULL;
-
+
return &solver_data->eps;
}
switch (state)
{
case ::dealii::SolverControl::iterate:
- solver_data->reason = EPS_CONVERGED_ITERATING;
- break;
-
+ solver_data->reason = EPS_CONVERGED_ITERATING;
+ break;
+
case ::dealii::SolverControl::success:
- solver_data->reason = static_cast<EPSConvergedReason>(1);
- break;
-
+ solver_data->reason = static_cast<EPSConvergedReason>(1);
+ break;
+
case ::dealii::SolverControl::failure:
- if (solver_control.last_step() > solver_control.max_steps())
- solver_data->reason = EPS_DIVERGED_ITS;
- else
- solver_data->reason = EPS_DIVERGED_BREAKDOWN;
- break;
-
+ if (solver_control.last_step() > solver_control.max_steps())
+ solver_data->reason = EPS_DIVERGED_ITS;
+ else
+ solver_data->reason = EPS_DIVERGED_BREAKDOWN;
+ break;
+
default:
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
}
}
int
SolverBase::convergence_test (EPS /*eps */,
- PetscScalar /*kr */,
- PetscScalar /*ki */,
- PetscReal /*residual_norm */,
- PetscReal */*estimated_error */,
+ PetscScalar /*kr */,
+ PetscScalar /*ki */,
+ PetscReal /*residual_norm */,
+ PetscReal */*estimated_error */,
void */*solver_control_x*/)
{
// This function is undefined (future reference only).
/* ---------------- Generalized Davidson ----------------- */
SolverGeneralizedDavidson::SolverGeneralizedDavidson (SolverControl &cn,
- const MPI_Comm &mpi_communicator,
- const AdditionalData &data)
+ const MPI_Comm &mpi_communicator,
+ const AdditionalData &data)
:
SolverBase (cn, mpi_communicator),
additional_data (data)
/* ------------------ Jacobi Davidson -------------------- */
SolverJacobiDavidson::SolverJacobiDavidson (SolverControl &cn,
- const MPI_Comm &mpi_communicator,
- const AdditionalData &data)
+ const MPI_Comm &mpi_communicator,
+ const AdditionalData &data)
:
SolverBase (cn, mpi_communicator),
additional_data (data)
dealii_sparse_matrix.get_sparsity_pattern();
if (matrix.get() == 0 ||
- m() != n_rows ||
+ m() != n_rows ||
n_nonzero_elements() != sparsity_pattern.n_nonzero_elements())
- {
- SparsityPattern trilinos_sparsity;
- trilinos_sparsity.reinit (input_row_map, input_col_map, sparsity_pattern);
- reinit (trilinos_sparsity);
- }
+ {
+ SparsityPattern trilinos_sparsity;
+ trilinos_sparsity.reinit (input_row_map, input_col_map, sparsity_pattern);
+ reinit (trilinos_sparsity);
+ }
// fill the values. the same as above: go through all rows of the matrix,
// and then all columns. since the sparsity patterns of the input matrix
// optimized diagonal
AssertDimension(it->column(), row);
if (std::fabs(it->value()) > drop_tolerance)
- {
- values[col] = it->value();
- row_indices[col++] = it->column();
- }
+ {
+ values[col] = it->value();
+ row_indices[col++] = it->column();
+ }
++select_index;
++it;
}
std::map<unsigned int, double> boundary_values;
for (std::set<unsigned int>::const_iterator p=boundary_dofs.begin();
- p != boundary_dofs.end(); ++p)
+ p != boundary_dofs.end(); ++p)
boundary_values[*p] = 0;
Vector<number> dummy(matrix.m());
std::map<unsigned int, double> boundary_values;
for (std::set<unsigned int>::const_iterator p=boundary_dofs.begin();
- p != boundary_dofs.end(); ++p)
+ p != boundary_dofs.end(); ++p)
boundary_values[*p] = 0;
BlockVector<number> dummy(matrix.n_block_rows());
// set entries of this line to zero except for the diagonal
// entry
- for (typename SparseMatrix<number>::iterator
- p = matrix.begin(dof_number);
- p != matrix.end(dof_number); ++p)
- if (p->column() != dof_number)
- p->value() = 0.;
+ for (typename SparseMatrix<number>::iterator
+ p = matrix.begin(dof_number);
+ p != matrix.end(dof_number); ++p)
+ if (p->column() != dof_number)
+ p->value() = 0.;
// set right hand side to
// wanted value: if main diagonal
// since that is the diagonal element and thus the present
// row
for (typename SparseMatrix<number>::iterator
- q = matrix.begin(dof_number)+1;
- q != matrix.end(dof_number); ++q)
+ q = matrix.begin(dof_number)+1;
+ q != matrix.end(dof_number); ++q)
{
const unsigned int row = q->column();
// (row,dof_number)
bool (*comp)(const typename SparseMatrix<number>::iterator::value_type p,
const unsigned int column)
- = &column_less_than<typename SparseMatrix<number>::iterator>;
+ = &column_less_than<typename SparseMatrix<number>::iterator>;
const typename SparseMatrix<number>::iterator
p = Utilities::lower_bound(matrix.begin(row)+1,
matrix.end(row),
// the diagonal block
for (typename SparseMatrix<number>::iterator
q = (block_index.first == block_row ?
- transpose_matrix.begin(block_index.second)+1 :
- transpose_matrix.begin(block_index.second));
- q != transpose_matrix.end(block_index.second);
- ++q)
+ transpose_matrix.begin(block_index.second)+1 :
+ transpose_matrix.begin(block_index.second));
+ q != transpose_matrix.end(block_index.second);
+ ++q)
{
// get the number
// of the column in
// sub-matrices
bool (*comp)(const typename SparseMatrix<number>::iterator::value_type p,
const unsigned int column)
- = &column_less_than<typename SparseMatrix<number>::iterator>;
+ = &column_less_than<typename SparseMatrix<number>::iterator>;
typename SparseMatrix<number>::iterator p = this_matrix.end();