Hello Rajat,
Thanks for the suggestion. It worked for the simple case of moving all
nodes of the triangulation by a constant displacement. In the actual
application, I want to move the mesh using the values from a
parallel:distributed displacement field which also has periodic boundary
conditions. For example
for (int d=0; d<3; ++d)
vertexDisplacement[d]=diisplacementVector(cell->vertex_dof_index(vertex_no,d));
However, this leads to illegal memory access on the displacement vector
when iterating over the ghost cells of the triangulation.
I have also attached below a minimal non-working example of my issue. It
would be very helpful if I could use a function on the lines of
triangulation.communicate_locally_moved_vertices(..).
Thank you,
Sambit
On Monday, December 4, 2017 at 2:22:03 PM UTC-5, RAJAT ARORA wrote:
>
> Hello Sambit,
>
> Can you try doing this?
> Also move the vertices of the ghost cell and avoid calling
> dftPtr->triangulation.communicate_locally_moved_vertices(locally_owned_vertices);
>
> When I tried to use
> triangulation.communicate_locally_moved_vertices(locally_owned_vertices)
> last time, something wierd (Dont exactly remember what, Code was not doing
> what I was expecting) happened with me. But Since then, I have avoided
> using it.
>
>
>
> Thanks.
>
> On Monday, December 4, 2017 at 1:22:37 PM UTC-5, Sambit Das wrote:
>>
>> Hello Dr. Arndt,
>>
>> Thank you for your reply. My apologies for not being clear on the "
>> breaks the periodic face pairs match". Following is the error message I get
>> when I run on parallel.
>>
>> *An error occurred in line <3699> of file
>> </home/vikramg/DFT-FE-softwares/softwareCentos/dealii/dealii-8.5.1/source/grid/grid_tools.cc>
>>
>> in function*
>> * void
>> dealii::GridTools::match_periodic_face_pairs(std::set<std::pair<CellIterator,
>>
>> unsigned int>, std::less<std::pair<CellIterator, unsigned int>>,
>> std::allocator<std::pair<CellIterator, unsigned int>>> &,
>> std::set<std::pair<dealii::identity<CellIterator>::type, unsigned int>,
>> std::less<std::pair<dealii::identity<CellIterator>::type, unsigned int>>,
>> std::allocator<std::pair<dealii::identity<CellIterator>::type, unsigned
>> int>>> &, int,
>> std::vector<dealii::GridTools::PeriodicFacePair<CellIterator>,
>> std::allocator<dealii::GridTools::PeriodicFacePair<CellIterator>>> &, const
>> dealii::Tensor<1, CellIterator::AccessorType::space_dimension, double> &,
>> const dealii::FullMatrix<double> &) [with CellIterator =
>> dealii::TriaIterator<dealii::CellAccessor<3, 3>>]*
>> *The violated condition was: *
>> * n_matches == pairs1.size() && pairs2.size() == 0*
>> *Additional information: *
>> * Unmatched faces on periodic boundaries*
>>
>>
>> I suspect this has something to do with the ghost nodes across the
>> periodic boundary not being handled correctly. I am right now creating a
>> minimal working example of my bug. I will post that soon.
>>
>> Best,
>> Sambit
>>
>>
>> On Monday, December 4, 2017 at 8:53:51 AM UTC-5, Daniel Arndt wrote:
>>>
>>> Sambit,
>>>
>>> I am trying to move all parallel triangulation nodes by a constant
>>>> displacement, but that breaks the periodic face pairs match when I
>>>> call GridTools::collect_periodic_faces(...). I use the following code for
>>>> the mesh movement. The dftPtr->triangulation has periodicity constraints
>>>> using add_periodicity(...).
>>>>
>>> What exactly do you mean by "that breaks the periodic face pairs match"?
>>> What is the error you are observing?
>>> Can you provide us with a minimal example that shows the problem so we
>>> can check?
>>>
>>> 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.
//
////C++ headers
////
#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/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <fstream>
#include <cmath>
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);
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);
std::vector<GridTools::PeriodicFacePair<typename parallel::distributed::Triangulation<3>::cell_iterator> > periodicity_vector2;
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_vector2);
}
std::cout << "number of elements: "
<< triangulation.n_global_active_cells()
<< std::endl;
}
