Hi, I was trying to change my code written with dealii 9.0.1 to make it work with dealii 9.2.0 where I am facing some problem.
One small issue which I noticed is the way how dealii constructs the
Constraint Matrix has probably changed. In order to illustrate this thing,
I have prepared an example (code along with input attached to the mail for
reference ) where I read a triangulation, Mark PBC Boundary IDs, add
periodicity, mark periodic faces and then print the constraint matrix info.
For code if compiled with dealii 9.0.1. the output is :
Periodic Constraints info:
0 8: 1
1 9: 1
4 10: 1
5 11: 1
12 16: 1
13 17: 1
Periodic Constraints info:
0 12: 1
1 13: 1
2 14: 1
3 15: 1
8 16: 1
9 17: 1
whereas same code compiled with dealii 9.2.0 and ran with same example gave
the following output:
Periodic Constraints info:
8 0: 1
9 1: 1
10 4: 1
11 5: 1
16 12: 1
17 13: 1
Periodic Constraints info:
12 0: 1
13 1: 1
14 2: 1
15 3: 1
16 8: 1
17 9: 1
The order in which Periodic info is received for plus and and minus faces
has reversed, May I know if this behavior was desired or is just a bug ?
Looking forward to some response.
Regards,
Aditya
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/d466a606-1248-473b-9798-53ccf9d3580an%40googlegroups.com.
#include <deal.II/distributed/tria.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <iostream>
#include <fstream>
using NumberType = double;
using UnsignedIntType = unsigned int;
template <UnsignedIntType dim>
using PeriodicPair1 = dealii::GridTools::PeriodicFacePair<typename
dealii::parallel::distributed::Triangulation<dim>::cell_iterator>;
template <UnsignedIntType dim>
using PeriodicPair2 = dealii::GridTools::PeriodicFacePair<typename
dealii::DoFHandler<dim>::cell_iterator>;
void Read_meshFile(const std::string meshFileName,
dealii::parallel::distributed::Triangulation<2>&
triangulation)
{
typename dealii::GridIn<2>::Format meshInputFormat =
dealii::GridIn<2>::abaqus;
dealii::GridIn<2> gridIn;
gridIn.attach_triangulation(triangulation);
std::ifstream gridInputStream(meshFileName);
gridIn.read(gridInputStream, meshInputFormat);
std::cout << "nElements: " << triangulation.n_active_cells()
<< "\t nNodes: " << triangulation.n_vertices()
<< std::endl;
}
void Mark_periodicFaces(
const dealii::types::boundary_id &boundaryIDStart,
const dealii::DoFHandler<2> & doFHandler,
std::array<std::vector<PeriodicPair2<2>>,2> & periodicFacePairs)
{
for (unsigned int boundaryFacePairCtr = 0;
boundaryFacePairCtr < 2;
++boundaryFacePairCtr)
{
// \note here plus faces should come first as later the constraints
// will be added for the first boundary_ID among the arguments
dealii::GridTools::collect_periodic_faces(
doFHandler,
100 + boundaryIDStart + boundaryFacePairCtr,
boundaryIDStart + boundaryFacePairCtr,
boundaryFacePairCtr,
periodicFacePairs[boundaryFacePairCtr]);
}
}
void Add_periodicityToTriangulation(
const dealii::types::boundary_id & boundaryIDStart,
dealii::parallel::distributed::Triangulation<2> &triangulation)
{
std::array<std::vector<PeriodicPair1<2>>,2> periodicFacePairs;
for (unsigned int boundaryFacePairCtr = 0;
boundaryFacePairCtr < 2;
++boundaryFacePairCtr)
{
dealii::GridTools::collect_periodic_faces(
triangulation,
100 + boundaryIDStart + boundaryFacePairCtr,
boundaryIDStart + boundaryFacePairCtr,
boundaryFacePairCtr,
periodicFacePairs[boundaryFacePairCtr]);
triangulation.add_periodicity(periodicFacePairs[boundaryFacePairCtr]);
}
}
void Mark_boundaryIDsPBC(dealii::parallel::distributed::Triangulation<2>
&triangulation, NumberType &rveEdgeLength,
NumberType &rveBoundaryCorrectionFactor,
dealii::types::boundary_id & pbcBoundaryIDStart)
{
const auto nFacesPerCell = dealii::GeometryInfo<2>::faces_per_cell;
unsigned countBoundaryFaces = 0, nBoundaryFacesExpected = 0;
NumberType maxFaceCenterCoord =
0.5 * rveEdgeLength;
const auto &halfRVELength =
0.5 * rveEdgeLength;
const auto &spatialToleranceCorrection =
rveBoundaryCorrectionFactor * rveEdgeLength;
for (const auto &cell : triangulation.active_cell_iterators())
{
if (cell->at_boundary())
{
for (unsigned int faceCtr = 0;
faceCtr < nFacesPerCell;
++faceCtr)
{
const auto ¤tFace = cell->face(faceCtr);
if (currentFace->at_boundary())
{
if (cell->is_locally_owned())
++nBoundaryFacesExpected;
for (dealii::types::boundary_id i = 0;
i < static_cast<dealii::types::boundary_id>(2);
++i)
{
const auto ¤tFaceCenter = currentFace->center();
// store most deviated boundary coordinate
if (std::fabs(currentFaceCenter[i]) >
maxFaceCenterCoord)
maxFaceCenterCoord =
std::fabs(currentFaceCenter[i]);
// mark face IDs
if (std::fabs(currentFaceCenter[i] + halfRVELength) <
(1e-8 +
spatialToleranceCorrection))
{
currentFace->set_boundary_id(
pbcBoundaryIDStart +
i);
if (cell->is_locally_owned())
++countBoundaryFaces;
}
if (std::fabs(currentFaceCenter[i] - halfRVELength) <
(1e-8 +
spatialToleranceCorrection))
{
currentFace->set_boundary_id(
100 +
pbcBoundaryIDStart +
i);
if (cell->is_locally_owned())
++countBoundaryFaces;
}
}
}
}
}
}
// sum all marked locally owned faces and raise assertion if the sum does
// not equal the expected number of boundary faces
countBoundaryFaces =
dealii::Utilities::MPI::sum(countBoundaryFaces, MPI_COMM_WORLD);
nBoundaryFacesExpected =
dealii::Utilities::MPI::sum(nBoundaryFacesExpected, MPI_COMM_WORLD);
double globalMaxFaceCenterCoord;
// get largest deviating coordinate to print in assertion
Assert(countBoundaryFaces == nBoundaryFacesExpected,
dealii::ExcMessage(
"Not all boundary faces were marked! Increase the boundary "
"correction factor! The maximum boundary face coordinate is " +
std::to_string(globalMaxFaceCenterCoord) + "."));
}
void Print_information(
const std::array<std::vector<PeriodicPair2<2>>,2> &
periodicBoundaryFacePairs)
{
// \note ComponentMask will be crucial for multifield formulations
dealii::ComponentMask componentMask;
dealii::ConstraintMatrix periodicConstraints;
for (unsigned int boundaryFacePairCtr = 0;
boundaryFacePairCtr < 2;
boundaryFacePairCtr++)
{
periodicConstraints.clear();
dealii::DoFTools::make_periodicity_constraints<dealii::DoFHandler<2>>(
periodicBoundaryFacePairs[boundaryFacePairCtr]
, periodicConstraints
,componentMask);
std::cout << std::endl << "Periodic Constraints info: " <<std::endl;
periodicConstraints.print(std::cout);
}
}
//main function
//
int main(int argc, char *argv[])
{
Assert(argc == 3 || argc == 4, dealii::ExcMessage("Insufficient number of
input arguments. "
"The program runs as \n ./verify-mesh meshFileName rveEdgeLength"));
dealii::Utilities::MPI::MPI_InitFinalize mpiInitialization(argc, argv, 1);
MPI_Comm const & mpiCommunicator(MPI_COMM_WORLD);
const UnsignedIntType dim = 2;
const std::string inputMeshFile = argv[1];
NumberType rveEdgeLength = std::stoi(argv[2]);
NumberType rveBoundaryCorrectionFactor;
dealii::types::boundary_id pbcBoundaryIDstart = 101;
if (argc == 4) {
rveBoundaryCorrectionFactor = std::stod(argv[3]);
AssertThrow((rveBoundaryCorrectionFactor <= 0.05) &&
(rveBoundaryCorrectionFactor >= 0.0),
dealii::ExcMessage("RVE boundary correction factor out of
range. Must be within [0.0, 0.05]"));
} else{
rveBoundaryCorrectionFactor = 0.015;
}
std::array<std::vector<PeriodicPair2<2>>,2> periodicFacePairs;
dealii::parallel::distributed::Triangulation<dim>
triangulation(mpiCommunicator, dealii::Triangulation<dim>::maximum_smoothing);
Read_meshFile(inputMeshFile, triangulation);
const dealii::FE_Q<dim, dim> fe_Q(1);
dealii::FESystem< dim, dim > feSystem(fe_Q, 2);
dealii::DoFHandler<dim,dim> dofHandler(triangulation);
dofHandler.distribute_dofs(feSystem);
Mark_boundaryIDsPBC(triangulation,rveEdgeLength,rveBoundaryCorrectionFactor,pbcBoundaryIDstart);
Add_periodicityToTriangulation(pbcBoundaryIDstart, triangulation);
Mark_periodicFaces(pbcBoundaryIDstart, dofHandler, periodicFacePairs);
Print_information(periodicFacePairs);
}
2D-testmodel.inp
Description: chemical/gamess-input
