diff --git a/include/dft.h b/include/dft.h index 884529008..9dfa0aa1a 100644 --- a/include/dft.h +++ b/include/dft.h @@ -229,7 +229,7 @@ namespace dftfe { * @brief generate image charges and update k point cartesian coordinates based * on current lattice vectors */ - void initImageChargesUpdateKPoints(); + void initImageChargesUpdateKPoints(bool flag=true); /** */ @@ -585,6 +585,8 @@ namespace dftfe { meshMovementGaussianClass d_gaussianMovePar; std::vector> d_gaussianMovementAtomsNetDisplacements; + std::vector > d_controlPointLocationsCurrentMove; + double d_gaussianConstantAutoMove; /// volume of the domain double d_domainVolume; diff --git a/include/meshMovement.h b/include/meshMovement.h index aaed85d30..e8d3a02a2 100644 --- a/include/meshMovement.h +++ b/include/meshMovement.h @@ -92,6 +92,17 @@ namespace dftfe { const std::vector > & controlPointDisplacements, const double controllingParameter, const bool moveSubdivided = false)=0; + + + + /*virtual std::pair moveMeshTwoStep(const std::vector > & controlPointLocations1, + const std::vector > & controlPointLocations2, + const std::vector > & controlPointDisplacements1, + const std::vector > & controlPointDisplacements2, + const double controllingParameter1, + const double controllingParameter2, + const bool moveSubdivided = false) = 0;*/ + virtual void computeIncrement()=0; /// vector of displacements of the triangulation vertices diff --git a/include/meshMovementGaussian.h b/include/meshMovementGaussian.h index 11777a418..97c3fbe31 100644 --- a/include/meshMovementGaussian.h +++ b/include/meshMovementGaussian.h @@ -50,11 +50,29 @@ namespace dftfe { const std::vector > & controlPointDisplacements, const double controllingParameter, const bool moveSubdivided = false); + + + + std::pair moveMeshTwoStep(const std::vector > & controlPointLocations1, + const std::vector > & controlPointLocations2, + const std::vector > & controlPointDisplacements1, + const std::vector > & controlPointDisplacements2, + const double controllingParameter1, + const double controllingParameter2, + const bool moveSubdivided = false); + private: /** @brief internal function which computes the nodal increment field in the local processor * */ void computeIncrement(); + + void computeIncrementTwoStep(const std::vector > & controlPointLocations1, + const std::vector > & controlPointLocations2, + const std::vector > & controlPointDisplacements1, + const std::vector > & controlPointDisplacements2, + const double controllingParameter1, + const double controllingParameter2); /// internal: storage for coordinates of the control points to which the Gaussians are attached std::vector > d_controlPointLocations; diff --git a/src/dft/dft.cc b/src/dft/dft.cc index f28daacd3..2aca6f4e6 100644 --- a/src/dft/dft.cc +++ b/src/dft/dft.cc @@ -447,7 +447,7 @@ namespace dftfe { // generate image charges and update k point cartesian coordinates based on current lattice vectors template - void dftClass::initImageChargesUpdateKPoints() + void dftClass::initImageChargesUpdateKPoints(bool flag) { TimerOutput::Scope scope (computing_timer, "image charges and k point generation"); pcout<<"-----------Simulation Domain bounding vectors (lattice vectors in fully periodic case)-------------"< periodicBc(3,false); periodicBc[0]=dftParameters::periodicX;periodicBc[1]=dftParameters::periodicY;periodicBc[2]=dftParameters::periodicZ; const double tol=1e-6; + + if(flag) + { for(unsigned int i = 0; i < atomLocationsFractional.size(); ++i) { for(unsigned int idim = 0; idim < 3; ++idim) @@ -484,6 +487,7 @@ namespace dftfe { "corresponding algorithm.")); } } + } generateImageCharges(d_pspCutOff, d_imageIds, diff --git a/src/dft/moveAtoms.cc b/src/dft/moveAtoms.cc index 2cc329a55..3075be0b2 100644 --- a/src/dft/moveAtoms.cc +++ b/src/dft/moveAtoms.cc @@ -13,7 +13,7 @@ // // --------------------------------------------------------------------- // -// @author Sambit Das +// @author Sambit Das and Phani Motamarri // namespace internal{ @@ -87,8 +87,8 @@ template void dftClass::updateAtomPositionsAndMoveMesh(const std::vector > & globalAtomsDisplacements) { const int numberGlobalAtoms = atomLocations.size(); - const int numberImageCharges = d_imageIds.size(); - const int totalNumberAtoms = numberGlobalAtoms + numberImageCharges; + int numberImageCharges = d_imageIds.size(); + int totalNumberAtoms = numberGlobalAtoms + numberImageCharges; std::vector latticeVectorsFlattened(9,0.0); for (unsigned int idim=0; idim<3; idim++) @@ -104,22 +104,26 @@ void dftClass::updateAtomPositionsAndMoveMesh(const std::vector periodicBc(3,false); periodicBc[0]=dftParameters::periodicX;periodicBc[1]=dftParameters::periodicY;periodicBc[2]=dftParameters::periodicZ; - std::vector > controlPointLocations; - std::vector > controlPointDisplacements; + std::vector > controlPointLocationsInitialMove; + std::vector > controlPointDisplacementsInitialMove; - std::vector > tempDispClosestTriaVerticesToAtoms; + std::vector > controlPointLocationsCurrentMove; + std::vector > controlPointDisplacementsCurrentMove; + + d_gaussianMovementAtomsNetDisplacements.resize(numberGlobalAtoms); + std::vector > tempGaussianMovementAtomsNetDisplacements; double maxDispAtom=-1; for(unsigned int iAtom=0;iAtom < numberGlobalAtoms; iAtom++) { - d_dispClosestTriaVerticesToAtoms[iAtom]+=globalAtomsDisplacements[iAtom]; - const double netDisp = d_dispClosestTriaVerticesToAtoms[iAtom].norm(); + d_gaussianMovementAtomsNetDisplacements[iAtom]+=globalAtomsDisplacements[iAtom]; + const double netDisp = d_gaussianMovementAtomsNetDisplacements[iAtom].norm(); if(netDisp>maxDispAtom) maxDispAtom=netDisp; } - tempDispClosestTriaVerticesToAtoms = d_dispClosestTriaVerticesToAtoms; + tempGaussianMovementAtomsNetDisplacements = d_gaussianMovementAtomsNetDisplacements; unsigned int useGaussian = 0; const double tol=1e-6; @@ -135,9 +139,13 @@ void dftClass::updateAtomPositionsAndMoveMesh(const std::vector atomCoor; int atomId=iAtom; @@ -190,7 +198,7 @@ void dftClass::updateAtomPositionsAndMoveMesh(const std::vector::updateAtomPositionsAndMoveMesh(const std::vector temp; int atomId; + Point<3> atomCoor; if(iAtom < numberGlobalAtoms) { - atomId = iAtom; - temp = d_closestTriaVertexToAtomsLocation[atomId]; - d_dispClosestTriaVerticesToAtoms.push_back(tempDispClosestTriaVerticesToAtoms[iAtom]); + d_gaussianMovementAtomsNetDisplacements.push_back(tempGaussianMovementAtomsNetDisplacements[iAtom]); } else { - Point<3> imageCoor; - Point<3> correspondingAtomCoor; - unsigned int iImage = iAtom - numberGlobalAtoms; - imageCoor[0] = d_imagePositions[iImage][0]; - imageCoor[1] = d_imagePositions[iImage][1]; - imageCoor[2] = d_imagePositions[iImage][2]; - const int atomId=d_imageIds[iImage]; - correspondingAtomCoor[0] = atomLocations[atomId][2]; - correspondingAtomCoor[1] = atomLocations[atomId][3]; - correspondingAtomCoor[2] = atomLocations[atomId][4]; - temp = d_closestTriaVertexToAtomsLocation[atomId]+(imageCoor-correspondingAtomCoor); - d_dispClosestTriaVerticesToAtoms.push_back(d_dispClosestTriaVerticesToAtoms[atomId]); + const int atomId=d_imageIds[iAtom-numberGlobalAtoms]; + d_gaussianMovementAtomsNetDisplacements.push_back(d_gaussianMovementAtomsNetDisplacements[atomId]); } - controlPointLocations.push_back(temp); - controlPointDisplacements.push_back(d_dispClosestTriaVerticesToAtoms[iAtom]); + controlPointLocationsInitialMove.push_back(d_closestTriaVertexToAtomsLocation[iAtom]); + controlPointDisplacementsInitialMove.push_back(d_dispClosestTriaVerticesToAtoms[iAtom]); + controlPointDisplacementsCurrentMove.push_back(d_gaussianMovementAtomsNetDisplacements[iAtom]); } @@ -300,9 +300,9 @@ void dftClass::updateAtomPositionsAndMoveMesh(const std::vectorgetGaussianGeneratorParameter()<<" as max displacement magnitude: "<< maxDispAtom<< " is below " << break1 <<" Bohr"< meshQualityMetrics=gaussianMove.moveMesh(controlPointLocations,controlPointDisplacements,forcePtr->getGaussianGeneratorParameter()); + const std::pair meshQualityMetrics=gaussianMove.moveMeshTwoStep(controlPointLocationsInitialMove,d_controlPointLocationsCurrentMove,controlPointDisplacementsInitialMove,controlPointDisplacementsCurrentMove,d_gaussianConstantAutoMove,forcePtr->getGaussianGeneratorParameter()); unsigned int autoMesh=0; - if (meshQualityMetrics.first || meshQualityMetrics.second>1.2*d_autoMeshMaxJacobianRatio) + if (meshQualityMetrics.first || meshQualityMetrics.second>1.3*d_autoMeshMaxJacobianRatio) autoMesh=1; MPI_Bcast(&(autoMesh), 1, @@ -314,7 +314,40 @@ void dftClass::updateAtomPositionsAndMoveMesh(const std::vectorgetGaussianGeneratorParameter()<getGaussianGeneratorParameter()<getGaussianGeneratorParameter()< atomCoor; + int atomId=iAtom; + atomCoor[0] = atomLocations[iAtom][2]; + atomCoor[1] = atomLocations[iAtom][3]; + atomCoor[2] = atomLocations[iAtom][4]; + + Point newCoord; + for(unsigned int idim=0; idim newFracCoord=internal::wrapAtomsAcrossPeriodicBc(newCoord, + corner, + latticeVectorsFlattened, + periodicBc); + //for synchrozination + MPI_Bcast(&(newFracCoord[0]), + 3, + MPI_DOUBLE, + 0, + MPI_COMM_WORLD); + + atomLocationsFractional[iAtom][2]=newFracCoord[0]; + atomLocationsFractional[iAtom][3]=newFracCoord[1]; + atomLocationsFractional[iAtom][4]=newFracCoord[2]; + } + + } + init(0); //for (unsigned int iAtom=0;iAtom ::moveMeshToAtoms(Triangulation<3,3> & triangulationMove, d_closestTriaVertexToAtomsLocation = closestTriaVertexToAtomsLocation; d_dispClosestTriaVerticesToAtoms = dispClosestTriaVerticesToAtoms; + d_controlPointLocationsCurrentMove.clear(); + for (unsigned int iAtom=0;iAtom atomCoor; + if(iAtom < numberGlobalAtoms) + { + atomCoor[0] = atomLocations[iAtom][2]; + atomCoor[1] = atomLocations[iAtom][3]; + atomCoor[2] = atomLocations[iAtom][4]; + } + else + { + atomCoor[0] = d_imagePositions[iAtom-numberGlobalAtoms][0]; + atomCoor[1] = d_imagePositions[iAtom-numberGlobalAtoms][1]; + atomCoor[2] = d_imagePositions[iAtom-numberGlobalAtoms][2]; + } + d_controlPointLocationsCurrentMove.push_back(atomCoor); + } + double minDist=1e+6; for (unsigned int i=0;i ::moveMeshToAtoms(Triangulation<3,3> & triangulationMove, dispClosestTriaVerticesToAtoms, gaussianConstant, moveSubdivided); + + d_gaussianConstantAutoMove = gaussianConstant; + timer_movemesh.exit_section("move mesh to atoms: move mesh"); AssertThrow(!meshQualityMetrics.first,ExcMessage("Negative jacobian created after moving closest nodes to atoms. Suggestion: increase refinement near atoms")); diff --git a/src/triangulation/meshMovement/meshMovementAffineTransform.cc b/src/triangulation/meshMovement/meshMovementAffineTransform.cc index 226b78f2b..a032de7f9 100644 --- a/src/triangulation/meshMovement/meshMovementAffineTransform.cc +++ b/src/triangulation/meshMovement/meshMovementAffineTransform.cc @@ -54,6 +54,7 @@ std::pair meshMovementAffineTransform::moveMesh(const std::vector

::vertices_per_cell; diff --git a/src/triangulation/meshMovement/meshMovementGaussian.cc b/src/triangulation/meshMovement/meshMovementGaussian.cc index 00f9f6e7b..69d4a51ba 100644 --- a/src/triangulation/meshMovement/meshMovementGaussian.cc +++ b/src/triangulation/meshMovement/meshMovementGaussian.cc @@ -54,6 +54,168 @@ std::pair meshMovementGaussianClass::moveMesh(const std::vector meshMovementGaussianClass::moveMeshTwoStep(const std::vector > & controlPointLocationsInitialMove, + const std::vector > & controlPointLocationsCurrentMove, + const std::vector > & controlPointDisplacementsInitialMove, + const std::vector > & controlPointDisplacementsCurrentMove, + const double controllingParameterInitialMove, + const double controllingParameterCurrentMove, + const bool moveSubdivided) +{ + //d_controlPointLocations = controlPointLocations; + //d_controlPointDisplacements=controlPointDisplacements; + //d_controllingParameter=controllingParameter; + //writeMesh("meshUnmoved.vtu"); + MPI_Barrier(mpi_communicator); + if (dftParameters::verbosity==2) + pcout << "Computing triangulation displacement increment caused by gaussian generator displacements..." << std::endl; + + initIncrementField(); + computeIncrementTwoStep(controlPointLocationsInitialMove, + controlPointLocationsCurrentMove, + controlPointDisplacementsInitialMove, + controlPointDisplacementsCurrentMove, + controllingParameterInitialMove, + controllingParameterCurrentMove); + finalizeIncrementField(); + if (dftParameters::verbosity==2) + pcout << "...Computed triangulation displacement increment" << std::endl; + if(moveSubdivided) + moveSubdividedMesh(); + + updateTriangulationVertices(); + std::pair returnData = movedMeshCheck(); + //writeMesh("meshMoved.vtu"); + return returnData; +} + + +void meshMovementGaussianClass::computeIncrementTwoStep(const std::vector > & controlPointLocationsInitialMove, + const std::vector > & controlPointLocationsCurrentMove, + const std::vector > & controlPointDisplacementsInitialMove, + const std::vector > & controlPointDisplacementsCurrentMove, + const double controllingParameterInitialMove, + const double controllingParameterCurrentMove) +{ + unsigned int vertices_per_cell = GeometryInfo::vertices_per_cell; + std::vector vertex_touched(d_dofHandlerMoveMesh.get_triangulation().n_vertices(), + false); + + std::vector > nodalCoordinatesUpdated(d_dofHandlerMoveMesh.get_triangulation().n_vertices()); + DoFHandler<3>::active_cell_iterator cell = d_dofHandlerMoveMesh.begin_active(), endc = d_dofHandlerMoveMesh.end(); + + for(; cell!=endc; ++cell) + if(!cell->is_artificial()) + for(unsigned int i=0; ivertex_index(i); + + if(vertex_touched[global_vertex_no]) + continue; + + Point nodalCoor = cell->vertex(i); + + if(!vertex_touched[global_vertex_no]) + nodalCoordinatesUpdated[global_vertex_no] = nodalCoor; + + vertex_touched[global_vertex_no]=true; + + + int overlappedControlPointId=-1; + + //check for case where control point locations coincide with nodal vertex locations + for(unsigned int jControl = 0; jControl < controlPointLocationsInitialMove.size(); jControl++) + { + const double distance=(nodalCoor-controlPointLocationsInitialMove[jControl]).norm(); + if (distance < 1e-5) + { + overlappedControlPointId=jControl; + break; + } + } + + for(unsigned int iControl=0;iControl vertex_dof_index(i,idim); + + if(!d_constraintsMoveMesh.is_constrained(globalDofIndex)) + { + d_incrementalDisplacement[globalDofIndex] + +=gaussianWeight*controlPointDisplacementsInitialMove[iControl][idim]; + + nodalCoordinatesUpdated[global_vertex_no][idim] += d_incrementalDisplacement[globalDofIndex]; + } + + } + + } + + } + + + DoFHandler<3>::active_cell_iterator cellStep2 = d_dofHandlerMoveMesh.begin_active(), endcStep2 = d_dofHandlerMoveMesh.end(); + std::vector vertex_touchedNew(d_dofHandlerMoveMesh.get_triangulation().n_vertices(), + false); + + for(; cellStep2!=endcStep2; ++cellStep2) + if(!cellStep2->is_artificial()) + for(unsigned int i = 0; i < vertices_per_cell; ++i) + { + const unsigned global_vertex_no = cellStep2->vertex_index(i); + + if(vertex_touchedNew[global_vertex_no]) + continue; + + vertex_touchedNew[global_vertex_no] = true; + int overlappedControlPointId=-1; + for(unsigned int jControl = 0; jControl < controlPointLocationsCurrentMove.size(); jControl++) + { + const double distance=(nodalCoordinatesUpdated[global_vertex_no]-controlPointLocationsCurrentMove[jControl]).norm(); + if (distance < 1e-5) + { + overlappedControlPointId=jControl; + break; + } + } + + for(unsigned int iControl = 0;iControl < controlPointLocationsCurrentMove.size(); iControl++) + { + if(overlappedControlPointId != iControl && overlappedControlPointId != -1) + continue; + + const double rsq=(nodalCoordinatesUpdated[global_vertex_no]-controlPointLocationsCurrentMove[iControl]).norm_square(); + const double gaussianWeight=dftParameters::reproducible_output? + std::exp(-(rsq)/std::pow(controllingParameterCurrentMove,2)) + :std::exp(-(rsq*rsq)/std::pow(controllingParameterCurrentMove,4)); + + for (unsigned int idim=0; idim < C_DIM ; idim++) + { + const unsigned int globalDofIndex = cellStep2->vertex_dof_index(i,idim); + + if(!d_constraintsMoveMesh.is_constrained(globalDofIndex)) + { + d_incrementalDisplacement[globalDofIndex] + +=gaussianWeight*controlPointDisplacementsCurrentMove[iControl][idim]; + + } + + } + + } + + } +} + //The triangulation nodes corresponding to control point location are constrained to only //their corresponding controlPointDisplacements. In other words for those nodes we don't consider overlapping @@ -78,36 +240,37 @@ void meshMovementGaussianClass::computeIncrement() Point nodalCoor = cell->vertex(i); int overlappedControlPointId=-1; - for (unsigned int jControl=0;jControl vertex_dof_index(i,idim); + { + const unsigned int globalDofIndex=cell->vertex_dof_index(i,idim); - if(!d_constraintsMoveMesh.is_constrained(globalDofIndex)) - d_incrementalDisplacement[globalDofIndex] - +=gaussianWeight*d_controlPointDisplacements[iControl][idim]; + if(!d_constraintsMoveMesh.is_constrained(globalDofIndex)) + d_incrementalDisplacement[globalDofIndex] + +=gaussianWeight*d_controlPointDisplacements[iControl][idim]; - } - } + } + } } }