This was caught by cppcheck.
number2 sum = 0.;
const size_type n_rows = m();
const number *val_ptr = &this->values[0];
- const number2 *v_ptr;
for (size_type row=0; row<n_rows; ++row)
{
number s = 0.;
const number *const val_end_of_row = val_ptr+n_rows;
- v_ptr = v.begin();
+ const number2 *v_ptr = v.begin();
while (val_ptr != val_end_of_row)
s += number(*val_ptr++) * number(*v_ptr++);
const size_type n_rows = m();
const size_type n_cols = n();
const number *val_ptr = &this->values[0];
- const number2 *v_ptr;
for (size_type row=0; row<n_rows; ++row)
{
number s = 0.;
const number *const val_end_of_row = val_ptr+n_cols;
- v_ptr = v.begin();
+ const number2 *v_ptr = v.begin();
while (val_ptr != val_end_of_row)
s += number(*val_ptr++) * number(*v_ptr++);
/* reinit *this to 0 */
this->reinit(A.m(), A.n());
- double SLik2 = 0.0, SLikLjk = 0.0;
for (size_type i=0; i< this->n_cols(); i++)
{
- SLik2 = 0.0;
+ double SLik2 = 0.0;
for (size_type j = 0; j < i; j++)
{
- SLikLjk = 0.0;
+ double SLikLjk = 0.0;
for (size_type k =0; k<j; k++)
{
SLikLjk += (*this)(i,k)*(*this)(j,k);
| additional_data.compute_condition_number
| additional_data.compute_all_condition_numbers
| additional_data.compute_eigenvalues;
- double eigen_beta_alpha = 0;
// vectors used for eigenvalue
// computations
try
{
+ double eigen_beta_alpha = 0;
+
// define some aliases for simpler access
VectorType &g = *Vr;
VectorType &d = *Vz;
d.reinit(x, true);
h.reinit(x, true);
- double gh,alpha,beta;
+ double gh,beta;
// compute residual. if vector is
// zero, then short-circuit the
it++;
A.vmult(h,d);
- alpha = d*h;
+ double alpha = d*h;
Assert(alpha != 0., ExcDivideByZero());
alpha = gh/alpha;
double r_l2 = 0;
double r0 = 0;
- double tau = 0;
+ double tau = 0;
double c = 0;
- double gamma = 0;
- double s = 0;
- double d_ = 0;
- double d = 0;
+ double s = 0;
+ double d_ = 0;
// The iteration step.
unsigned int j = 1;
A.vmult(*u[2],v);
u[2]->add (-std::sqrt(delta[1]/delta[0]), *u[0]);
- gamma = *u[2] * v;
+ const double gamma = *u[2] * v;
u[2]->add (-gamma / std::sqrt(delta[1]), *u[1]);
*m[0] = v;
e[1] = -c * std::sqrt(delta[2]);
}
- d = std::sqrt (d_*d_ + delta[2]);
+ const double d = std::sqrt (d_*d_ + delta[2]);
if (j>1)
tau *= s / c;
int it=0;
- double tau, rho, theta=0, sigma, alpha, psi, theta_old, rho_old, beta;
+ double tau, rho, theta=0;
double res;
d.reinit(x);
// Step 1
A.vmult(t,q);
// Step 2
- sigma = q*t;
+ const double sigma = q*t;
//TODO:[?] Find a really good breakdown criterion. The absolute one detects breakdown instead of convergence
if (std::fabs(sigma/rho) < additional_data.breakdown)
return IterationResult(SolverControl::iterate, std::fabs(sigma/rho));
// Step 3
- alpha = rho/sigma;
+ const double alpha = rho/sigma;
v.add(-alpha,t);
// Step 4
- theta_old = theta;
+ const double theta_old = theta;
theta = v*v/tau;
- psi = 1./(1.+theta);
+ const double psi = 1./(1.+theta);
tau *= theta*psi;
d.sadd(psi*theta_old, psi*alpha, p);
if (std::fabs(rho) < additional_data.breakdown)
return IterationResult(SolverControl::iterate, std::fabs(rho));
// Step 7
- rho_old = rho;
+ const double rho_old = rho;
precondition.vmult(q,v);
rho = q*v;
- beta = rho/rho_old;
+ const double beta = rho/rho_old;
p.sadd(beta,v);
precondition.vmult(q,p);
}
const unsigned int n_owned = (vector_partitioner->local_range().second-
vector_partitioner->local_range().first);
const std::size_t n_ghosts = ghost_dofs.size();
- unsigned int n_unique_ghosts= 0;
#ifdef DEBUG
for (std::vector<unsigned int>::iterator dof = dof_indices.begin();
dof!=dof_indices.end(); ++dof)
IndexSet ghost_indices (vector_partitioner->size());
if (n_ghosts > 0)
{
+ unsigned int n_unique_ghosts = 0;
// since we need to go back to the local_to_global indices and
// replace the temporary numbering of ghosts by the real number in
// the index set, we need to store these values
std::vector<internal::color::PartitionWork<Worker>*> worker(n_workers);
std::vector<internal::color::PartitionWork<Worker>*> blocked_worker(n_blocked_workers);
unsigned int worker_index = 0, slice_index = 0;
- unsigned int spawn_index = 0, spawn_index_new = 0;
+ unsigned int spawn_index = 0;
int spawn_index_child = -2;
internal::MPIComCompress<OutVector> *worker_compr = new(root->allocate_child())
internal::MPIComCompress<OutVector>(dst);
for (unsigned int part=0;
part<task_info.partition_color_blocks_row_index.size()-1; part++)
{
- spawn_index_new = worker_index;
+ const unsigned int spawn_index_new = worker_index;
if (part == 0)
worker[worker_index] = new(worker_compr->allocate_child())
internal::color::PartitionWork<Worker>(func,slice_index,task_info,false);
// all vector data, then handle the
// scalar data sets that have been
// left over
- unsigned int i, pt_data_vector_dim;
+ unsigned int i;
std::string vector_name;
for (i=0; i<data_filter.n_data_sets(); ++i)
{
// Allocate space for the point data
// Must be either 1D or 3D
- pt_data_vector_dim = data_filter.get_data_set_dim(i);
+ const unsigned int pt_data_vector_dim = data_filter.get_data_set_dim(i);
vector_name = data_filter.get_data_set_name(i);
// Create the dataspace for the point data
(void)component;
Assert (component==0, ExcIndexRange(component,0,1));
- double prod;
double sum = 0;
for (unsigned int monom = 0; monom < exponents.n_rows(); ++monom)
{
- prod = 1;
+ double prod = 1;
for (unsigned int s=0; s< dim; ++s)
{
if (p[s] < 0)
{
const double d = this->center.distance(points[k]);
const double r = this->radius;
- double e = 0.;
double val = 0.;
if (d<this->radius)
{
- e = -r*r/(r*r-d*d);
+ const double e = -r*r/(r*r-d*d);
if (e>-50)
val = numbers::E * exp(e);
}
long double z = std::cos(numbers::PI * (i-.25)/(n+.5));
long double pp;
- long double p1, p2, p3;
+ long double p1;
// Newton iteration
do
{
// compute L_n (z)
p1 = 1.;
- p2 = 0.;
+ long double p2 = 0.;
for (unsigned int j=0; j<n; ++j)
{
- p3 = p2;
+ const long double p3 = p2;
p2 = p1;
p1 = ((2.*j+1.)*z*p2-j*p3)/(j+1);
}
for (unsigned int i=0; i<m; ++i)
x[i] = - std::cos( (long double) (2*i+1)/(2*m) * numbers::PI );
- long double r, s, J_x, f, delta;
+ long double s, J_x, f, delta;
for (unsigned int k=0; k<m; ++k)
{
- r = x[k];
+ long double r = x[k];
if (k>0)
r = (r + x[k-1])/2;
// the Jacobi polynomial is evaluated
// using a recursion formula.
std::vector<long double> p(n+1);
- int v, a1, a2, a3, a4;
// initial values P_0(x), P_1(x):
p[0] = 1.0L;
for (unsigned int i=1; i<=(n-1); ++i)
{
- v = 2*i + alpha + beta;
- a1 = 2*(i+1)*(i + alpha + beta + 1)*v;
- a2 = (v + 1)*(alpha*alpha - beta*beta);
- a3 = v*(v + 1)*(v + 2);
- a4 = 2*(i+alpha)*(i+beta)*(v + 2);
+ const int v = 2*i + alpha + beta;
+ const int a1 = 2*(i+1)*(i + alpha + beta + 1)*v;
+ const int a2 = (v + 1)*(alpha*alpha - beta*beta);
+ const int a3 = v*(v + 1)*(v + 2);
+ const int a4 = 2*(i+alpha)*(i+beta)*(v + 2);
p[i+1] = static_cast<long double>( (a2 + a3*x)*p[i] - a4*p[i-1])/a1;
} // for
double eps = 1e-8;
unsigned int q_id = 0; // Current quad point index.
- double area = 0;
Tensor<1,2> dist;
for (unsigned int box=0; box<4; ++box)
{
dist = (singularity-GeometryInfo<2>::unit_cell_vertex(box));
dist = Point<2>(std::abs(dist[0]), std::abs(dist[1]));
- area = dist[0]*dist[1];
+ double area = dist[0]*dist[1];
if (area > eps)
for (unsigned int q=0; q<quads[box].size(); ++q, ++q_id)
{
typename dealii::internal::p4est::types<dim>::quadrant &quad)
{
int i, l = quad.level;
- int child_id;
dealii::types::global_dof_index dealii_index =
triangulation->get_p4est_tree_to_coarse_cell_permutation()[treeidx];
for (i = 0; i < l; i++)
{
typename dealii::Triangulation<dim,spacedim>::cell_iterator cell (triangulation, i, dealii_index);
- child_id = dealii::internal::p4est::functions<dim>::quadrant_ancestor_id (&quad, i + 1);
+ const int child_id = dealii::internal::p4est::functions<dim>::quadrant_ancestor_id (&quad, i + 1);
Assert (cell->has_children (), ExcMessage ("p4est quadrant does not correspond to a cell!"));
dealii_index = cell->child_index(child_id);
}
}
int factorial_i = 1;
- int factorial_ij;
- int factorial_j;
for (int i = 2; i < (int) this->dofs_per_face; ++i)
{
interpolation_matrix (i, i) = std::pow (0.5, i);
factorial_i *= i;
- factorial_j = factorial_i;
- factorial_ij = 1;
+ int factorial_j = factorial_i;
+ int factorial_ij = 1;
for (int j = i + 1; j < (int) this->dofs_per_face; ++j)
{
interpolation_matrix (1, 1) = 1.0;
int factorial_i = 1;
- int factorial_ij;
- int factorial_j;
for (int i = 2; i < (int) this->dofs_per_face; ++i)
{
interpolation_matrix (i, i) = std::pow (0.5, i);
factorial_i *= i;
- factorial_j = factorial_i;
- factorial_ij = 1;
+ int factorial_j = factorial_i;
+ int factorial_ij = 1;
for (int j = i + 1; j < (int) this->dofs_per_face; ++j)
{
interpolation_matrix (3, vertex) = 0.25;
int factorial_i = 1;
- int factorial_ij;
- int factorial_j;
- int factorial_k;
- int factorial_kl;
- int factorial_l;
for (int i = 2; i <= (int) this->degree; ++i)
{
j = j + 2;
}
- factorial_k = 1;
+ int factorial_k = 1;
for (int j = 2; j <= (int) this->degree; ++j)
{
interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j);
factorial_k *= j;
- factorial_kl = 1;
- factorial_l = factorial_k;
+ int factorial_kl = 1;
+ int factorial_l = factorial_k;
for (int k = j + 1; k < (int) this->degree; ++k)
{
}
factorial_i *= i;
- factorial_j = factorial_i;
- factorial_ij = 1;
+ int factorial_j = factorial_i;
+ int factorial_ij = 1;
for (int j = i + 1; j <= (int) this->degree; ++j)
{
{
interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k);
factorial_k *= k;
- factorial_l = factorial_k;
- factorial_kl = 1;
+ int factorial_l = factorial_k;
+ int factorial_kl = 1;
for (int l = k + 1; l <= (int) this->degree; ++l)
{
{
interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k);
factorial_k *= k;
- factorial_l = factorial_k;
- factorial_kl = 1;
+ int factorial_l = factorial_k;
+ int factorial_kl = 1;
for (int l = k + 1; l <= (int) this->degree; ++l)
{
interpolation_matrix (2, vertex) = 0.25;
int factorial_i = 1;
- int factorial_ij;
- int factorial_j;
- int factorial_k;
- int factorial_kl;
- int factorial_l;
for (int i = 2; i <= (int) this->degree; ++i)
{
interpolation_matrix (i + source_fe.degree + 1, i + this->degree + 1) = tmp;
interpolation_matrix (i + 2 * source_fe.degree, i + 2 * this->degree) = tmp;
factorial_i *= i;
- factorial_j = factorial_i;
- factorial_ij = 1;
+ int factorial_j = factorial_i;
+ int factorial_ij = 1;
for (int j = i + 1; j <= (int) this->degree; ++j)
{
factorial_j *= j;
tmp = std::pow (0.5, j) * factorial_j / (factorial_i * factorial_ij);
interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp;
- factorial_k = 1;
+ int factorial_k = 1;
for (int k = 2; k <= (int) this->degree; ++k)
{
interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k);
factorial_k *= k;
- factorial_l = factorial_k;
- factorial_kl = 1;
+ int factorial_l = factorial_k;
+ int factorial_kl = 1;
for (int l = k + 1; l <= (int) this->degree; ++l)
{
}
}
- factorial_k = 1;
+ int factorial_k = 1;
for (int j = 2; j <= (int) this->degree; ++j)
{
interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j);
factorial_k *= j;
- factorial_l = factorial_k;
- factorial_kl = 1;
+ int factorial_l = factorial_k;
+ int factorial_kl = 1;
for (int k = j + 1; k <= (int) this->degree; ++k)
{
interpolation_matrix (1, vertex) = 0.25;
int factorial_i = 1;
- int factorial_ij;
- int factorial_j;
- int factorial_k;
- int factorial_kl;
- int factorial_l;
for (int i = 2; i <= (int) this->degree; ++i)
{
j = j + 2;
}
- factorial_k = 1;
+ int factorial_k = 1;
for (int j = 2; j <= (int) this->degree; ++j)
{
interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j);
factorial_k *= j;
- factorial_l = factorial_k;
- factorial_kl = 1;
+ int factorial_kl = 1;
+ int factorial_l = factorial_k;
for (int k = j + 1; k <= (int) this->degree; ++k)
{
}
factorial_i *= i;
- factorial_j = factorial_i;
- factorial_ij = 1;
+ int factorial_j = factorial_i;
+ int factorial_ij = 1;
for (int j = i + 1; j <= (int) this->degree; ++j)
{
interpolation_matrix (i + 2 * source_fe.degree, j + 3 * this->degree - 1) = tmp;
tmp *= 2.0;
interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp;
- factorial_k = 1;
+ int factorial_k = 1;
for (int k = 2; k <= (int) this->degree; ++k)
{
interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k);
factorial_k *= k;
- factorial_l = factorial_k;
- factorial_kl = 1;
+ int factorial_l = factorial_k;
+ int factorial_kl = 1;
for (int l = k + 1; l <= (int) this->degree; ++l)
{
interpolation_matrix (3, 3) = 1.0;
int factorial_i = 1;
- int factorial_ij;
- int factorial_j;
- int factorial_k;
- int factorial_kl;
- int factorial_l;
for (int i = 2; i <= (int) this->degree; ++i)
{
tmp *= -1.0;
interpolation_matrix (i + source_fe.degree + 1, i + this->degree + 1) = tmp;
interpolation_matrix (i + 3 * source_fe.degree - 1, i + 3 * this->degree - 1) = tmp;
- factorial_k = 1;
+ int factorial_k = 1;
for (int j = 2; j <= (int) this->degree; ++j)
{
interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j);
factorial_k *= j;
- factorial_l = factorial_k;
- factorial_kl = 1;
+ int factorial_l = factorial_k;
+ int factorial_kl = 1;
for (int k = j + 1; k <= (int) this->degree; ++k)
{
}
factorial_i *= i;
- factorial_j = factorial_i;
- factorial_ij = 1;
+ int factorial_j = factorial_i;
+ int factorial_ij = 1;
for (int j = i + 1; j <= (int) this->degree; ++j)
{
tmp *= 2.0;
interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = tmp;
interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp;
- factorial_k = 1;
+ int factorial_k = 1;
for (int k = 2; k <= (int) this->degree; ++k)
{
interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k);
factorial_k *= k;
- factorial_l = factorial_k;
- factorial_kl = 1;
+ int factorial_l = factorial_k;
+ int factorial_kl = 1;
for (int l = k + 1; l <= (int) this->degree; ++l)
{
// Surface can be created from ELSET, or directly from cells
// If elsets_list contains a key with specific name - refers to that ELSET, otherwise refers to cell
std::istringstream iss (line);
- char comma;
int el_idx;
int face_number;
char temp;
else
{
// Surface refers directly to elements
+ char comma;
iss >> el_idx >> comma >> temp >> face_number;
quad_node_list = get_global_node_numbers (el_idx, face_number);
quad_node_list.insert (quad_node_list.begin(), b_indicator);
// draw the legend
if (svg_flags.draw_legend) out << '\n' << " <!-- legend -->" << '\n';
- unsigned int line_offset = 0;
-
additional_width = 0;
if (!svg_flags.margin) additional_width = static_cast<unsigned int>(.5 + (height/100.) * 2.5);
// explanation of the cell labeling
if (svg_flags.draw_legend && (svg_flags.label_level_number || svg_flags.label_cell_index || svg_flags.label_material_id || svg_flags.label_subdomain_id || svg_flags.label_level_subdomain_id ))
{
+ unsigned int line_offset = 0;
out << " <rect x=\"" << width + additional_width << "\" y=\"" << static_cast<unsigned int>(.5 + (height/100.) * margin_in_percent)
<< "\" width=\"" << static_cast<unsigned int>(.5 + (height/100.) * (40. - margin_in_percent)) << "\" height=\"" << static_cast<unsigned int>(.5 + height * .165) << "\"/>" << '\n';
typename std::vector<typename Container::cell_iterator>::const_iterator uniform_cell;
for (uniform_cell=uniform_cells.begin(); uniform_cell!=uniform_cells.end(); ++uniform_cell)
{
- bool repeat_vertex;
for (unsigned int v=0; v<GeometryInfo<Container::dimension>::vertices_per_cell; ++v)
{
Point<Container::space_dimension> position=(*uniform_cell)->vertex (v);
- repeat_vertex=false;
+ bool repeat_vertex=false;
for (unsigned int m=0; m<i; ++m)
{
const std::pair<size_type, size_type>
local_range = rhs.local_range();
- int ierr;
// Try to copy all the rows of the matrix one by one. In case of error
// (i.e., the column indices are different), we need to abort and blow
// away the matrix.
int n_entries, rhs_n_entries;
TrilinosScalar *value_ptr, *rhs_value_ptr;
int *index_ptr, *rhs_index_ptr;
- ierr = rhs.matrix->ExtractMyRowView (row_local, rhs_n_entries,
- rhs_value_ptr, rhs_index_ptr);
+ int ierr = rhs.matrix->ExtractMyRowView (row_local, rhs_n_entries,
+ rhs_value_ptr, rhs_index_ptr);
(void)ierr;
Assert (ierr == 0, ExcTrilinosError(ierr));
local_range = rhs.local_range();
const bool same_col_map = matrix->ColMap().SameAs(rhs.matrix->ColMap());
- int ierr;
for (size_type row=local_range.first; row < local_range.second; ++row)
{
const int row_local =
int n_entries, rhs_n_entries;
TrilinosScalar *value_ptr, *rhs_value_ptr;
int *index_ptr, *rhs_index_ptr;
- ierr = rhs.matrix->ExtractMyRowView (row_local, rhs_n_entries,
- rhs_value_ptr, rhs_index_ptr);
+ int ierr = rhs.matrix->ExtractMyRowView (row_local, rhs_n_entries,
+ rhs_value_ptr, rhs_index_ptr);
(void)ierr;
Assert (ierr == 0, ExcTrilinosError(ierr));
// otherwise first flush Trilinos caches
sparsity_pattern->compress ();
- // get a representation of the present row
- int ncols;
-
colnum_cache.reset (new std::vector<size_type> (sparsity_pattern->row_length(this->a_row)));
if (colnum_cache->size() > 0)
{
- int ierr;
- ierr = sparsity_pattern->graph->ExtractGlobalRowCopy((TrilinosWrappers::types::int_type)this->a_row,
- colnum_cache->size(),
- ncols,
- (TrilinosWrappers::types::int_type *)&(*colnum_cache)[0]);
+ // get a representation of the present row
+ int ncols;
+ const int ierr = sparsity_pattern->graph->ExtractGlobalRowCopy
+ ((TrilinosWrappers::types::int_type)this->a_row,
+ colnum_cache->size(),
+ ncols,
+ (TrilinosWrappers::types::int_type *)&(*colnum_cache)[0]);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
AssertThrow (static_cast<std::vector<size_type>::size_type>(ncols) == colnum_cache->size(),
ExcInternalError());
Assert(Inters.IsDone(), ExcMessage("Could not project point."));
double minDistance = 1e7;
- double distance;
Point<3> result;
for (int i=0; i<Inters.NbPnt(); ++i)
{
- distance = point(origin).Distance(Inters.Pnt(i+1));
+ const double distance = point(origin).Distance(Inters.Pnt(i+1));
//cout<<"Point "<<i<<": "<<point(Inters.Pnt(i+1))<<" distance: "<<distance<<endl;
if (distance < minDistance)
{