test_propagators_00003.cpp
Example of multi-stepped vector field propagation on homogeneous hexa volume mesh with slip boundaries.Using: PropagateVectorField
To run: ./test_propagators_00003
visit: mimmo website
/*---------------------------------------------------------------------------*\
*
* mimmo
*
* Copyright (C) 2015-2021 OPTIMAD engineering Srl
*
* -------------------------------------------------------------------------
* License
* This file is part of mimmo.
*
* mimmo is free software: you can redistribute it and/or modify it
* under the terms of the GNU Lesser General Public License v3 (LGPL)
* as published by the Free Software Foundation.
*
* mimmo is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with mimmo. If not, see <http://www.gnu.org/licenses/>.
*
\ *---------------------------------------------------------------------------*/
#include "mimmo_propagators.hpp"
#if MIMMO_ENABLE_MPI
#include "mimmo_parallel.hpp"
#endif
// =================================================================================== //
// =================================================================================== //
// =================================================================================== //
mimmo::MimmoSharedPointer<mimmo::MimmoObject> createTestVolumeMesh(int rank, mimmo::MimmoSharedPointer<mimmo::MimmoObject> & boundary){
std::array<double,3> center({{0.0,0.0,0.0}});
double radiusin(2.0), radiusout(5.0);
double azimuthin(0.0), azimuthout(0.5*BITPIT_PI);
double heightbottom(-1.0), heighttop(1.0);
int nr(6), nt(10), nh(6);
double deltar = (radiusout - radiusin)/ double(nr);
double deltat = (azimuthout - azimuthin)/ double(nt);
double deltah = (heighttop - heightbottom)/ double(nh);
std::vector<std::array<double,3> > verts ((nr+1)*(nt+1)*(nh+1));
int counter = 0;
for(int k=0; k<=nh; ++k){
for(int j=0; j<=nt; ++j){
for(int i=0; i<=nr; ++i){
verts[counter][0] =(radiusin + i*deltar)*std::cos(azimuthin + j*deltat);
verts[counter][1] =(radiusin + i*deltar)*std::sin(azimuthin + j*deltat);
verts[counter][2] =(heightbottom + k*deltah);
++counter;
}
}
}
//create the volume mesh mimmo.
mimmo::MimmoSharedPointer<mimmo::MimmoObject> mesh = mimmo::MimmoSharedPointer<mimmo::MimmoObject>(new mimmo::MimmoObject(2));
if (rank == 0){
mesh->getPatch()->reserveVertices((nr+1)*(nt+1)*(nh+1));
//pump up the vertices
for(const auto & vertex : verts){
mesh->addVertex(vertex); //automatic id assigned to vertices.
}
mesh->getPatch()->reserveCells(nr*nt*nh);
//create connectivities for hexa elements
std::vector<long> conn(8,0);
for(int k=0; k<nh; ++k){
for(int j=0; j<nt; ++j){
for(int i=0; i<nr; ++i){
conn[0] = (nr+1)*(nt+1)*k + (nr+1)*j + i;
conn[1] = (nr+1)*(nt+1)*k + (nr+1)*j + i+1;
conn[2] = (nr+1)*(nt+1)*k + (nr+1)*(j+1) + i+1;
conn[3] = (nr+1)*(nt+1)*k + (nr+1)*(j+1) + i;
conn[4] = (nr+1)*(nt+1)*(k+1) + (nr+1)*j + i;
conn[5] = (nr+1)*(nt+1)*(k+1) + (nr+1)*j + i+1;
conn[6] = (nr+1)*(nt+1)*(k+1) + (nr+1)*(j+1) + i+1;
conn[7] = (nr+1)*(nt+1)*(k+1) + (nr+1)*(j+1) + i;
mesh->addConnectedCell(conn, bitpit::ElementType::HEXAHEDRON);
}
}
}
}// end rank 0
mesh->update();
boundary = mesh->extractBoundaryMesh();
//decompose manually the six boundary faces.
std::vector<std::vector<long>> bverts(6, std::vector<long>());
//vertex at azimuthal limits.
for(int k=0; k<=nh; ++k){
for(int i=0; i<=nr; ++i){
bverts[0].push_back((nr+1)*(nt+1)*k + i);
bverts[1].push_back((nr+1)*(nt+1)*k + (nr+1)*nt + i);
}
}
//vertex at quote limits.
for(int j=0; j<=nt; ++j){
for(int i=0; i<=nr; ++i){
bverts[2].push_back((nr+1)*j + i);
bverts[3].push_back((nr+1)*(nt+1)*nh + (nr+1)*j + i);
}
}
//vertex at radius limits.
for(int k=0; k<=nh; ++k){
for(int j=0; j<=nt; ++j){
bverts[4].push_back((nr+1)*(nt+1)*k + (nr+1)*j);
bverts[5].push_back((nr+1)*(nt+1)*k + (nr+1)*j + nr);
}
}
//put together
mimmo::MimmoPiercedVector<long> pidfaces;
long pidder = 1;
for(livector1D & vert : bverts ){
livector1D boundaryfaces = boundary->getCellFromVertexList(vert, true);
for(long id : boundaryfaces){
pidfaces.insert(id, pidder);
}
++pidder;
}
bverts.clear();
for(auto it= pidfaces.begin(); it!= pidfaces.end(); ++it){
boundary->setPIDCell(it.getId(), *it);
}
boundary->updateAdjacencies();
boundary->update();
boundary->getPatch()->write("piddedBoundary_test3");
return mesh;
}
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ //
mimmo::MimmoSharedPointer<mimmo::MimmoObject> pidExtractor(mimmo::MimmoSharedPointer<mimmo::MimmoObject> boundary, const livector1D & pids){
if(!boundary) return nullptr;
livector1D cellsID = boundary->extractPIDCells(pids);
livector1D vertID = boundary->getVertexFromCellList(cellsID);
extracted->getPatch()->reserveVertices(vertID.size());
extracted->getPatch()->reserveCells(cellsID.size());
for(long id: vertID){
extracted->addVertex(boundary->getVertexCoords(id), id);
}
for(long id: cellsID){
extracted->addCell(boundary->getCells().at(id), id);
}
extracted->updateAdjacencies();
extracted->update();
return extracted;
}
// =================================================================================== //
int test1() {
#if MIMMO_ENABLE_MPI
// Initialize mpi
int nProcs;
int rank;
MPI_Comm_size(MPI_COMM_WORLD, &nProcs);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
#else
int rank = 0;
#endif
mimmo::MimmoSharedPointer<mimmo::MimmoObject> mesh = createTestVolumeMesh(rank, bDirMesh);
#if MIMMO_ENABLE_MPI
/* Instantiation of a Partition object with default patition method space filling curve.
* Plot Optional results during execution active for Partition block.
*/
mimmo::Partition* partition = new mimmo::Partition();
partition->setPlotInExecution(true);
partition->setGeometry(mesh);
partition->setBoundaryGeometry(bDirMesh);
partition->setPartitionMethod(mimmo::PartitionMethod::PARTGEOM);
partition->exec();
#endif
mimmo::MimmoSharedPointer<mimmo::MimmoObject> bDIR = pidExtractor(bDirMesh, {{1,2}});
mimmo::MimmoSharedPointer<mimmo::MimmoObject> bSLIP = pidExtractor(bDirMesh, {{5,6}});
//TESTING THE VECTOR PROPAGATOR //////
// and the scalar field of Dirichlet values on pidded surface 0 and 1 nodes.
mimmo::MimmoPiercedVector<std::array<double,3>> bc_surf_3Dfield;
bc_surf_3Dfield.setGeometry(bDIR);
bc_surf_3Dfield.setDataLocation(mimmo::MPVLocation::POINT);
bc_surf_3Dfield.reserve(bDIR->getNVertices());
std::vector<long> d1 = bDIR->getVertexFromCellList(bDIR->extractPIDCells(1));
for(long id : d1){
bc_surf_3Dfield.insert(id, {{0.0, 0.0, 0.0}});
}
//apply a rotation of alpha around z centered in the origin 0,0,0, so that deformation is
//x-directed and equal to tan(alpha) * y-coord of the point.
//BEWARE IF THE ROTATION ARC IS similar to a linear segment, this deformation works,
//for large rotation, where the arc path is very different from a different segment, this
// could lead to not properly expected results.
double alpha = 25.0*BITPIT_PI/180.0;
std::vector<long> d2 = bDIR->getVertexFromCellList(bDIR->extractPIDCells(2));
for(long id : d2){
darray3E coord = bDIR->getVertexCoords(id);
double normc = std::sqrt(coord[0]*coord[0] + coord[1]*coord[1]);
bc_surf_3Dfield.insert(id, {{std::sin(alpha)*normc,
(std::cos(alpha)-1.0)*normc,
0.0}} );
}
// Now create a PropagateScalarField and solve the laplacian.
mimmo::PropagateVectorField * prop3D = new mimmo::PropagateVectorField();
prop3D->setName("test00003_PropagateVectorField");
prop3D->setGeometry(mesh);
prop3D->addDirichletBoundaryPatch(bDIR);
prop3D->addDirichletConditions(&bc_surf_3Dfield);
prop3D->addSlipBoundarySurface(bSLIP);
prop3D->addSlipReferenceSurface(bSLIP);
prop3D->setDamping(true);
prop3D->setDampingType(0);
prop3D->setDampingDecayFactor(2.0);
prop3D->setDampingInnerDistance(0.5);
prop3D->setDampingOuterDistance(2.5);
prop3D->setSolverMultiStep(40);
prop3D->setPlotInExecution(true);
prop3D->exec();
//deform the mesh and write deformation
mimmo::dmpvecarr3E * values3D = prop3D->getPropagatedField();
darray3E work;
for(auto it=values3D->begin(); it!=values3D->end(); ++it){
work = mesh->getVertexCoords(it.getId());
mesh->modifyVertex(work + *it, it.getId());
}
mesh->update();
mesh->getPatch()->write("test00003_deformedMesh");
bool check = false;
long targetNode = (10 +1)*(6+1)*3 + (6+1)*5 + 3;
darray3E value;
value.fill(-1.0*std::numeric_limits<double>::max());
if(values3D->exists(targetNode)) value = values3D->at(targetNode);
#if MIMMO_ENABLE_MPI
MPI_Allreduce(MPI_IN_PLACE, value.data(), 3, MPI_DOUBLE, MPI_MAX, prop3D->getCommunicator());
#endif
if(rank == 0){
std::cout<<targetNode<<" "<<value<<std::endl;
check = check || (norm2(values3D->at(targetNode)-std::array<double,3>({{0.601253, -0.284719, 0.}})) > 1.0E-5);
}
delete prop3D;
#if MIMMO_ENABLE_MPI
delete partition;
#endif
return check;
}
// =================================================================================== //
int main( int argc, char *argv[] ) {
BITPIT_UNUSED(argc);
BITPIT_UNUSED(argv);
#if MIMMO_ENABLE_MPI
MPI_Init(&argc, &argv);
#endif
int val = 1;
try{
val = test1() ;
}
catch(std::exception & e){
std::cout<<"test_propagators_00003 exited with an error of type : "<<e.what()<<std::endl;
return 1;
}
#if MIMMO_ENABLE_MPI
MPI_Finalize();
#endif
return val;
}
MimmoObject is the basic geometry container for mimmo library.
Definition: MimmoObject.hpp:143
void setGeometry(MimmoSharedPointer< MimmoObject > geo)
Definition: MimmoPiercedVector.tpp:314
@ POINT
1.8.17