Skip to content

Commit

Permalink
Merge pull request #233 from landinjm/range_based_loops
Browse files Browse the repository at this point in the history 
Switching to range based loops
  • Loading branch information
@landinjm
landinjm authored Sep 5, 2024
2 parents 22352a7 + fc66010 commit f66aa56
Show file tree
Hide file tree
Showing 17 changed files with 140 additions and 212 deletions.
4 changes: 1 addition & 3 deletions applications/CHiMaD_benchmark6b/customPDE.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -205,9 +205,7 @@ customPDE<dim, degree>::makeTriangulation(
}

// Mark the boundaries
typename parallel::distributed::Triangulation<dim>::cell_iterator cell4 = tria.begin(),
endc4 = tria.end();
for (; cell4 != endc4; ++cell4)
for (const auto &cell4 : tria.active_cell_iterators())
{
// Mark all of the faces
for (unsigned int face_number = 0; face_number < GeometryInfo<dim>::faces_per_cell;
Expand Down
31 changes: 14 additions & 17 deletions applications/_nucleating_precipitates_MgRE/equations.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -281,22 +281,19 @@ customPDE<dim, degree>::seedNucleus(
dealii::AlignedVector<dealii::VectorizedArray<double>> &source_terms,
dealii::VectorizedArray<double> &gamma) const
{
for (typename std::vector<nucleus<dim>>::const_iterator thisNucleus =
this->nuclei.begin();
thisNucleus != this->nuclei.end();
++thisNucleus)
for (const auto &thisNucleus : this->nuclei)
{
if (thisNucleus->seededTime + thisNucleus->seedingTime > this->currentTime)
if (thisNucleus.seededTime + thisNucleus.seedingTime > this->currentTime)
{
// Calculate the weighted distance function to the order parameter
// freeze boundary (weighted_dist = 1.0 on that boundary)
dealii::VectorizedArray<double> weighted_dist =
this->weightedDistanceFromNucleusCenter(
thisNucleus->center,
userInputs.get_nucleus_freeze_semiaxes(thisNucleus->orderParameterIndex),
userInputs.get_nucleus_rotation_matrix(thisNucleus->orderParameterIndex),
thisNucleus.center,
userInputs.get_nucleus_freeze_semiaxes(thisNucleus.orderParameterIndex),
userInputs.get_nucleus_rotation_matrix(thisNucleus.orderParameterIndex),
q_point_loc,
thisNucleus->orderParameterIndex);
thisNucleus.orderParameterIndex);

for (unsigned i = 0; i < gamma.size(); i++)
{
Expand All @@ -306,13 +303,13 @@ customPDE<dim, degree>::seedNucleus(

// Seed a nucleus if it was added to the list of nuclei this
// time step
if (thisNucleus->seedingTimestep == this->currentIncrement)
if (thisNucleus.seedingTimestep == this->currentIncrement)
{
// Set the rotation matrix for this nucleus (two possible
// per order parameter)
dealii::Tensor<2, dim, double> nucleus_rot_matrix;
nucleus_rot_matrix = userInputs.get_nucleus_rotation_matrix(
thisNucleus->orderParameterIndex);
thisNucleus.orderParameterIndex);

// Find the weighted distance to the outer edge of the
// nucleus and use it to calculate the order parameter
Expand All @@ -323,26 +320,26 @@ customPDE<dim, degree>::seedNucleus(
q_point_loc_element(j) = q_point_loc(j)[i];
}
double r = this->weightedDistanceFromNucleusCenter(
thisNucleus->center,
userInputs.get_nucleus_semiaxes(thisNucleus->orderParameterIndex),
thisNucleus.center,
userInputs.get_nucleus_semiaxes(thisNucleus.orderParameterIndex),
nucleus_rot_matrix,
q_point_loc_element,
thisNucleus->orderParameterIndex);
thisNucleus.orderParameterIndex);

double avg_semiaxis = 0.0;
for (unsigned int j = 0; j < dim; j++)
{
avg_semiaxis += thisNucleus->semiaxes[j];
avg_semiaxis += thisNucleus.semiaxes[j];
}
avg_semiaxis /= dim;

if (thisNucleus->orderParameterIndex == 1)
if (thisNucleus.orderParameterIndex == 1)
{
source_terms[0][i] =
0.5 *
(1.0 - std::tanh(avg_semiaxis * (r - 1.0) / interface_coeff));
}
else if (thisNucleus->orderParameterIndex == 2)
else if (thisNucleus.orderParameterIndex == 2)
{
source_terms[1][i] =
0.5 *
Expand Down
19 changes: 6 additions & 13 deletions applications/allenCahn_conserved/customPDE.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -242,14 +242,11 @@ customPDE<dim, degree>::solveIncrement(bool skip_time_dependent)
// applied, replace the ones that do with the Dirichlet value
// This clears the residual where we want to apply Dirichlet
// BCs, otherwise the solver sees a positive residual
for (std::map<types::global_dof_index, double>::const_iterator it =
this->valuesDirichletSet[fieldIndex]->begin();
it != this->valuesDirichletSet[fieldIndex]->end();
++it)
for (const auto &it : *this->valuesDirichletSet[fieldIndex])
{
if (this->residualSet[fieldIndex]->in_local_range(it->first))
if (this->residualSet[fieldIndex]->in_local_range(it.first))
{
(*this->residualSet[fieldIndex])(it->first) = 0.0;
(*this->residualSet[fieldIndex])(it.first) = 0.0;
}
}

Expand Down Expand Up @@ -341,16 +338,12 @@ customPDE<dim, degree>::solveIncrement(bool skip_time_dependent)

this->computeNonexplicitRHS();

for (std::map<types::global_dof_index,
double>::const_iterator it =
this->valuesDirichletSet[fieldIndex]->begin();
it != this->valuesDirichletSet[fieldIndex]->end();
++it)
for (const auto &it : *this->valuesDirichletSet[fieldIndex])
{
if (this->residualSet[fieldIndex]->in_local_range(
it->first))
it.first))
{
(*this->residualSet[fieldIndex])(it->first) = 0.0;
(*this->residualSet[fieldIndex])(it.first) = 0.0;
}
}

Expand Down
23 changes: 10 additions & 13 deletions applications/nucleationModel/equations.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,21 +133,18 @@ customPDE<dim, degree>::seedNucleus(
dealii::VectorizedArray<double> &source_term,
dealii::VectorizedArray<double> &gamma) const
{
for (typename std::vector<nucleus<dim>>::const_iterator thisNucleus =
this->nuclei.begin();
thisNucleus != this->nuclei.end();
++thisNucleus)
for (const auto &thisNucleus : this->nuclei)
{
if (thisNucleus->seededTime + thisNucleus->seedingTime > this->currentTime)
if (thisNucleus.seededTime + thisNucleus.seedingTime > this->currentTime)
{
// Calculate the weighted distance function to the order parameter
// freeze boundary (weighted_dist = 1.0 on that boundary)
dealii::VectorizedArray<double> weighted_dist =
this->weightedDistanceFromNucleusCenter(
thisNucleus->center,
userInputs.get_nucleus_freeze_semiaxes(thisNucleus->orderParameterIndex),
thisNucleus.center,
userInputs.get_nucleus_freeze_semiaxes(thisNucleus.orderParameterIndex),
q_point_loc,
thisNucleus->orderParameterIndex);
thisNucleus.orderParameterIndex);

for (unsigned i = 0; i < gamma.size(); i++)
{
Expand All @@ -157,7 +154,7 @@ customPDE<dim, degree>::seedNucleus(

// Seed a nucleus if it was added to the list of nuclei this
// time step
if (thisNucleus->seedingTimestep == this->currentIncrement)
if (thisNucleus.seedingTimestep == this->currentIncrement)
{
// Find the weighted distance to the outer edge of the
// nucleus and use it to calculate the order parameter
Expand All @@ -168,15 +165,15 @@ customPDE<dim, degree>::seedNucleus(
q_point_loc_element(j) = q_point_loc(j)[i];
}
double r = this->weightedDistanceFromNucleusCenter(
thisNucleus->center,
userInputs.get_nucleus_semiaxes(thisNucleus->orderParameterIndex),
thisNucleus.center,
userInputs.get_nucleus_semiaxes(thisNucleus.orderParameterIndex),
q_point_loc_element,
thisNucleus->orderParameterIndex);
thisNucleus.orderParameterIndex);

double avg_semiaxis = 0.0;
for (unsigned int j = 0; j < dim; j++)
{
avg_semiaxis += thisNucleus->semiaxes[j];
avg_semiaxis += thisNucleus.semiaxes[j];
}
avg_semiaxis /= dim;

Expand Down
23 changes: 10 additions & 13 deletions applications/nucleationModel_preferential/equations.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,21 +134,18 @@ customPDE<dim, degree>::seedNucleus(
dealii::VectorizedArray<double> &gamma) const
{
// Loop through all of the seeded nuclei
for (typename std::vector<nucleus<dim>>::const_iterator thisNucleus =
this->nuclei.begin();
thisNucleus != this->nuclei.end();
++thisNucleus)
for (const auto &thisNucleus : this->nuclei)
{
if (thisNucleus->seededTime + thisNucleus->seedingTime > this->currentTime)
if (thisNucleus.seededTime + thisNucleus.seedingTime > this->currentTime)
{
// Calculate the weighted distance function to the order parameter
// freeze boundary (weighted_dist = 1.0 on that boundary)
dealii::VectorizedArray<double> weighted_dist =
this->weightedDistanceFromNucleusCenter(
thisNucleus->center,
userInputs.get_nucleus_freeze_semiaxes(thisNucleus->orderParameterIndex),
thisNucleus.center,
userInputs.get_nucleus_freeze_semiaxes(thisNucleus.orderParameterIndex),
q_point_loc,
thisNucleus->orderParameterIndex);
thisNucleus.orderParameterIndex);

for (unsigned i = 0; i < gamma.size(); i++)
{
Expand All @@ -158,7 +155,7 @@ customPDE<dim, degree>::seedNucleus(

// Seed a nucleus if it was added to the list of nuclei this
// time step
if (thisNucleus->seedingTimestep == this->currentIncrement)
if (thisNucleus.seedingTimestep == this->currentIncrement)
{
// Find the weighted distance to the outer edge of the
// nucleus and use it to calculate the order parameter
Expand All @@ -169,15 +166,15 @@ customPDE<dim, degree>::seedNucleus(
q_point_loc_element(j) = q_point_loc(j)[i];
}
double r = this->weightedDistanceFromNucleusCenter(
thisNucleus->center,
userInputs.get_nucleus_semiaxes(thisNucleus->orderParameterIndex),
thisNucleus.center,
userInputs.get_nucleus_semiaxes(thisNucleus.orderParameterIndex),
q_point_loc_element,
thisNucleus->orderParameterIndex);
thisNucleus.orderParameterIndex);

double avg_semiaxis = 0.0;
for (unsigned int j = 0; j < dim; j++)
{
avg_semiaxis += thisNucleus->semiaxes[j];
avg_semiaxis += thisNucleus.semiaxes[j];
}
avg_semiaxis /= dim;

Expand Down
32 changes: 12 additions & 20 deletions src/OrderParameterRemapper/OrderParameterRemapper.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,17 +17,14 @@ OrderParameterRemapper<dim>::remap(
double transfer_buffer =
std::max(0.0, grain_representations.at(g).getDistanceToNeighbor() / 2.0);

typename dealii::DoFHandler<dim>::active_cell_iterator di =
dof_handler.begin_active();

// For now I have two loops, one where I copy the values from the old
// order parameter to the new one and a second where I zero out the
// old order parameter. This separation prevents writing zero-out
// values to the new order parameter. There probably is a more
// efficient way of doing this.
while (di != dof_handler.end())
for (const auto &dof : dof_handler.active_cell_iterators())
{
if (di->is_locally_owned())
if (dof->is_locally_owned())
{
unsigned int op_new = grain_representations.at(g).getOrderParameterId();
unsigned int op_old =
Expand All @@ -40,7 +37,7 @@ OrderParameterRemapper<dim>::remap(
v < dealii::GeometryInfo<dim>::vertices_per_cell;
v++)
{
if (di->vertex(v).distance(
if (dof->vertex(v).distance(
grain_representations.at(g).getCenter()) >
grain_representations.at(g).getRadius() + transfer_buffer)
{
Expand All @@ -56,22 +53,19 @@ OrderParameterRemapper<dim>::remap(
std::vector<dealii::types::global_dof_index> dof_indices(
dofs_per_cell,
0);
di->get_dof_indices(dof_indices);
dof->get_dof_indices(dof_indices);
for (unsigned int i = 0; i < dof_indices.size(); i++)
{
(*solution_fields.at(op_new))[dof_indices.at(i)] =
(*solution_fields.at(op_old))[dof_indices.at(i)];
}
}
}
++di;
}

di = dof_handler.begin_active();

while (di != dof_handler.end())
for (const auto &dof : dof_handler.active_cell_iterators())
{
if (di->is_locally_owned())
if (dof->is_locally_owned())
{
unsigned int op_new = grain_representations.at(g).getOrderParameterId();
unsigned int op_old =
Expand All @@ -84,7 +78,7 @@ OrderParameterRemapper<dim>::remap(
v < dealii::GeometryInfo<dim>::vertices_per_cell;
v++)
{
if (di->vertex(v).distance(
if (dof->vertex(v).distance(
grain_representations.at(g).getCenter()) >
grain_representations.at(g).getRadius() + transfer_buffer)
{
Expand All @@ -99,15 +93,14 @@ OrderParameterRemapper<dim>::remap(
std::vector<dealii::types::global_dof_index> dof_indices(
dofs_per_cell,
0);
di->get_dof_indices(dof_indices);
dof->get_dof_indices(dof_indices);

for (unsigned int i = 0; i < dof_indices.size(); i++)
{
(*solution_fields.at(op_old))[dof_indices.at(i)] = 0.0;
}
}
}
++di;
}
}
}
Expand Down Expand Up @@ -141,9 +134,9 @@ OrderParameterRemapper<dim>::remap_from_index_field(
// order parameter. This separation prevents writing zero-out values to
// the new order parameter. There probably is a more efficient way of
// doing this.
while (di != dof_handler.end())
for (const auto &dof : dof_handler.active_cell_iterators())
{
if (di->is_locally_owned())
if (dof->is_locally_owned())
{
unsigned int op_new = grain_representations.at(g).getOrderParameterId();
unsigned int op_old = grain_representations.at(g).getOldOrderParameterId();
Expand All @@ -153,7 +146,7 @@ OrderParameterRemapper<dim>::remap_from_index_field(
for (unsigned int v = 0; v < dealii::GeometryInfo<dim>::vertices_per_cell;
v++)
{
if (di->vertex(v).distance(grain_representations.at(g).getCenter()) >
if (dof->vertex(v).distance(grain_representations.at(g).getCenter()) >
grain_representations.at(g).getRadius() + transfer_buffer)
{
in_grain = false;
Expand All @@ -167,7 +160,7 @@ OrderParameterRemapper<dim>::remap_from_index_field(
{
std::vector<dealii::types::global_dof_index> dof_indices(dofs_per_cell,
0);
di->get_dof_indices(dof_indices);
dof->get_dof_indices(dof_indices);
for (unsigned int i = 0; i < dof_indices.size(); i++)
{
if (std::abs((*grain_index_field)[dof_indices.at(i)] -
Expand All @@ -179,7 +172,6 @@ OrderParameterRemapper<dim>::remap_from_index_field(
}
}
}
++di;
}
}
}
Expand Down
11 changes: 2 additions & 9 deletions src/matrixfree/boundaryConditions.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,10 +39,6 @@ MatrixFreePDE<dim, degree>::applyNeumannBCs()
if (userInputs.BC_list[starting_BC_list_index].var_BC_type[direction] ==
NEUMANN)
{
typename DoFHandler<dim>::active_cell_iterator cell = dofHandlersSet[0]
->begin_active(),
endc =
dofHandlersSet[0]->end();
FESystem<dim> *fe = FESet[currentFieldIndex];
QGaussLobatto<dim - 1> face_quadrature_formula(degree + 1);
FEFaceValues<dim> fe_face_values(*fe,
Expand All @@ -54,7 +50,7 @@ MatrixFreePDE<dim, degree>::applyNeumannBCs()
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

// Loop over each face on a boundary
for (; cell != endc; ++cell)
for (const auto &cell : dofHandlersSet[0]->active_cell_iterators())
{
for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
{
Expand Down Expand Up @@ -357,10 +353,7 @@ MatrixFreePDE<dim, degree>::setRigidBodyModeConstraints(
unsigned int vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;

// Loop over each locally owned cell
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler->begin_active(),
endc = dof_handler->end();

for (; cell != endc; ++cell)
for (const auto &cell : dof_handler->active_cell_iterators())
{
if (cell->is_locally_owned())
{
Expand Down
Loading

0 comments on commit f66aa56

Please sign in to comment.