Hello Dr. Arndt,
Thank you for your explanation. I followed your suggestions to create a
periodic match check without calling GridTools::collect_periodic_faces(..).
I use the PeriodicFacePairs
<https://www.google.com/url?q=https%3A%2F%2Fwww.dealii.org%2Fdeveloper%2Fdoxygen%2Fdeal.II%2FstructGridTools_1_1PeriodicFacePair.html&sa=D&sntz=1&usg=AFQjCNEz8WY6PTmmUdjh6DM9pZTIYQ-xGw>
called
before refinement or mesh movement, and query the periodic match with
orthogonal_equality(..). When I tested this I found that the new periodic
face match sanity check is satisfied on serial when I move the mesh, but on
parallel it fails after mesh movement i.e orthogonal_equality(..) returns
false for all three directions for a given PeriodicFacePair. I have
attached the minimal working example of the same (For the working example
to compile successfully the following patch needs to
used https://github.com/dealii/dealii/pull/5549/files )
This means I am still not moving the mesh consistently even after calling
"triangulation.communicate_locally_moved_vertices(locally_owned_vertices)".
Thank you,
Sambit
On Wednesday, December 6, 2017 at 3:37:49 PM UTC-5, Daniel Arndt wrote:
>
> Sambit,
>
> I am using the GridTools::collect_periodic_faces(..) as a sanity check
>> after moving the triangulation. I donot again set periodicity constraints.
>> The documentation also mentions "it is possible to call this function
>> several times with different boundary ids to generate a vector with all
>> periodic pairs". Moreover, in my actual application I never do refinement
>> as I read a pre-generated mesh file, where the similar error occurs after
>> moving the triangulation.
>>
> The problem with a distributed mesh is that in general not all processes
> know about the locations of all the vertices after the mesh has been moved.
> Luckily, it should be sufficient to call
> GridTools::collect_periodic_faces(..) only once for the coarse mesh before
> moving since only matching faces are stored.
> Calling add_periodicity
> <https://dealii.org/developer/doxygen/deal.II/classparallel_1_1distributed_1_1Triangulation.html#a7a3f626abb92aca7040eaa2c160686a0>
> beforehand
> should make sure that all the information on the locations of vertices and
> degrees of freedom across periodic boundaries is available to all relevant
> processes. This is important if you intend to use
> make_periodicity_constraints.
>
>
>> I did two more checks in the minimal example-
>> 1) by not calling GridTools::collect_periodic_faces(..) after
>> refinement, I do not get any error messages. But is there a way to check
>> whether the periodic match still holds in the moved triangulation without
>> calling GridTools::collect_periodic_faces(..)?
>>
> You can of course just inspect the content of the
> std::vector< PeriodicFacePair< typename MeshType::cell_iterator > > &
> matched_pairs,
> It stores PeriodicFacePairs
> <https://www.google.com/url?q=https%3A%2F%2Fwww.dealii.org%2Fdeveloper%2Fdoxygen%2Fdeal.II%2FstructGridTools_1_1PeriodicFacePair.html&sa=D&sntz=1&usg=AFQjCNEz8WY6PTmmUdjh6DM9pZTIYQ-xGw>
> that
> contain two cell iterator and the relevant face numbers. From these
> information you can obtain for example the center of the two faces.
> Since no positions are stored, matched_pairs should still be valid after
> moving the mesh if done consistently.
>
>
>> 2) by placing the GridTools::collect_periodic_faces(..) before moving the
>> triangulation but after refinement, it worked fine on serial and parallel,
>> which suggests something breaking after movement.
>>
> Yes, see above.
>
> Best,
> Daniel
>
--
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
//Include all deal.II header file
#include <deal.II/base/quadrature.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/table.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/base/tensor_function.h>
#include <deal.II/base/point.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/lapack_full_matrix.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_bicgstab.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/lac/parallel_vector.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/lac/slepc_solver.h>
#include <deal.II/base/config.h>
#include <deal.II/base/smartpointer.h>
//Include generic C++ headers
#include <fstream>
#include <iostream>
using namespace dealii;
int main (int argc, char *argv[])
{
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
const double L=20;
parallel::distributed::Triangulation<3> triangulation(MPI_COMM_WORLD);
GridGenerator::hyper_cube (triangulation, -L, L);
//mark faces
typename parallel::distributed::Triangulation<3>::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end();
for(; cell!=endc; ++cell)
{
for(unsigned int f=0; f < GeometryInfo<3>::faces_per_cell; ++f)
{
const Point<3> face_center = cell->face(f)->center();
if(cell->face(f)->at_boundary())
{
if (std::abs(face_center[0]+(L))<1.0e-8)
cell->face(f)->set_boundary_id(1);
else if (std::abs(face_center[0]-(L))<1.0e-8)
cell->face(f)->set_boundary_id(2);
else if (std::abs(face_center[1]+(L))<1.0e-8)
cell->face(f)->set_boundary_id(3);
else if (std::abs(face_center[1]-(L))<1.0e-8)
cell->face(f)->set_boundary_id(4);
else if (std::abs(face_center[2]+(L))<1.0e-8)
cell->face(f)->set_boundary_id(5);
else if (std::abs(face_center[2]-(L))<1.0e-8)
cell->face(f)->set_boundary_id(6);
}
}
}
std::vector<GridTools::PeriodicFacePair<typename parallel::distributed::Triangulation<3>::cell_iterator> > periodicity_vector;
for (int i = 0; i < 3; ++i)
{
GridTools::collect_periodic_faces(triangulation, /*b_id1*/ 2*i+1, /*b_id2*/ 2*i+2,/*direction*/ i, periodicity_vector);
}
triangulation.add_periodicity(periodicity_vector);
triangulation.refine_global(2);
for(unsigned int i=0; i< periodicity_vector.size(); ++i){{
std::vector<bool> isPeriodicFace(3);
for(unsigned int idim=0; idim<3; ++idim){
isPeriodicFace[idim]=GridTools::orthogonal_equality(periodicity_vector[i].cell[0]->face(periodicity_vector[i].face_idx[0]),periodicity_vector[i].cell[1]->face(periodicity_vector[i].face_idx[1]),idim);
}
if (isPeriodicFace[0]==false && isPeriodicFace[1]==false && isPeriodicFace[2]==false)
std::cout<<" Periodic Face pairs not matching periodically for any directions "<<std::endl;
} }
std::vector<bool> vertex_moved(triangulation.n_vertices(), false);
const std::vector<bool> locally_owned_vertices = GridTools::get_locally_owned_vertices(triangulation);
cell = triangulation.begin_active();
for(; cell!=endc; ++cell)
{
if (cell->is_locally_owned())
{
for (unsigned int vertex_no=0; vertex_no<GeometryInfo<3>::vertices_per_cell;
++vertex_no)
{
const unsigned global_vertex_no = cell->vertex_index(vertex_no);
if (vertex_moved[global_vertex_no]
|| !locally_owned_vertices[global_vertex_no])
continue;
Point<3> vertexDisplacement; vertexDisplacement[0]=1e-4; vertexDisplacement[1]=0; vertexDisplacement[2]=0;
cell->vertex(vertex_no) += vertexDisplacement;
vertex_moved[global_vertex_no] = true;
}
}
}
triangulation.communicate_locally_moved_vertices(locally_owned_vertices);
for(unsigned int i=0; i< periodicity_vector.size(); ++i)
{
std::vector<bool> isPeriodicFace(3);
for(unsigned int idim=0; idim<3; ++idim){
isPeriodicFace[idim]=GridTools::orthogonal_equality(periodicity_vector[i].cell[0]->face(periodicity_vector[i].face_idx[0]),periodicity_vector[i].cell[1]->face(periodicity_vector[i].face_idx[1]),idim);
}
if (isPeriodicFace[0]==false && isPeriodicFace[1]==false && isPeriodicFace[2]==false)
std::cout<<" Previously periodic matched face pairs not matching periodically for any directions after mesh movement "<<std::endl;
}
std::cout << "number of elements: "
<< triangulation.n_global_active_cells()
<< std::endl;
}
