37#include "PeriodicNodesOnlyMesh.hpp"
39template<
unsigned SPACE_DIM>
47 for (
unsigned i = 0; i < SPACE_DIM; i++)
56template<
unsigned SPACE_DIM>
63 EXCEPTION(
"Too many processors for the periodic domain width and cut off length. NumProcs should be less than or equal to Y (in 2D) or Z (in 3D) width / Cut off length.\n");
68 for (
unsigned i=0; i < SPACE_DIM; i++)
70 if (mIsDimPeriodic[i])
73 if ( fmod(mWidth[i]+1e-14,cutOffLength ) > 2e-14 )
75 EXCEPTION(
"The periodic width must be a multiple of cut off length.");
78 else if ( mWidth[i]/cutOffLength <= (2.0+1e-14) )
80 EXCEPTION(
"The periodic domain width cannot be less than 2*CutOffLength." );
84 domainSize[ i*2 ] = 0;
85 domainSize[ i*2+1 ] = mWidth[i];
90 this->AddNodesToBoxes();
93template<
unsigned SPACE_DIM>
99template<
unsigned SPACE_DIM>
103 assert(rDimension < SPACE_DIM);
105 if ( !mIsDimPeriodic[rDimension] )
112 width = mWidth[rDimension];
117template<
unsigned SPACE_DIM>
121 assert(!concreteMove);
123 c_vector<double, SPACE_DIM>& point_coord = point.
rGetLocation();
126 for (
unsigned i = 0; i < SPACE_DIM; i++ )
128 if ( mIsDimPeriodic[i] )
130 double local_point_coord = point_coord(i);
131 double local_width = mWidth[i];
134 double fudge_factor = 1.0e-13;
135 if ( local_point_coord >= local_width )
137 double new_coord = local_point_coord - local_width;
141 else if ( local_point_coord > (local_width-fudge_factor))
147 else if ( local_point_coord < 0.0 )
149 double new_coord = local_point_coord + local_width;
153 if ( new_coord > local_width-fudge_factor )
155 new_coord = local_width-fudge_factor;
163 this->GetNode(nodeIndex)->SetPoint(point);
166template<
unsigned SPACE_DIM>
174 SetNode(node_index, new_node_point);
180template<
unsigned SPACE_DIM>
183 c_vector<double, SPACE_DIM> vector = rLocation2 - rLocation1;
187 for (
unsigned i = 0; i < SPACE_DIM; i++ )
189 if ( mIsDimPeriodic[i] )
191 vector(i) = fmod(vector(i),mWidth[i]);
192 if ( vector(i) > 0.5*mWidth[i] )
194 vector(i) -= mWidth[i];
196 else if (vector(i) < -0.5*mWidth[i] )
198 vector(i) += mWidth[i];
207template<
unsigned SPACE_DIM>
212 for (
typename std::vector<
Node<SPACE_DIM> *>::iterator it_node = this->mNodes.begin();
213 it_node != this->mNodes.end(); ++it_node )
215 c_vector<double,SPACE_DIM> & location = (*it_node)->rGetModifiableLocation();
217 for (
unsigned i = 0; i < SPACE_DIM; i++ )
219 if ( mIsDimPeriodic[i] )
221 if ( location(i) < 0.0 )
223 location(i) = fmod( location(i), mWidth[i] ) + mWidth[i];
225 else if ( location(i) >= mWidth[i] )
227 location(i) = fmod( location(i), mWidth[i] );
#define EXCEPTION(message)
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
virtual void RefreshMesh()
virtual double GetWidth(const unsigned &rDimension) const
c_vector< double, DIM > & rGetLocation()
void SetCoordinate(unsigned i, double value)
ChastePoint< SPACE_DIM > GetPoint() const
unsigned AddNode(Node< SPACE_DIM > *pNewNode)
void SetUpBoxCollection(const std::vector< Node< SPACE_DIM > * > &rNodes)
c_vector< bool, SPACE_DIM > mIsDimPeriodic
void SetNode(unsigned nodeIndex, ChastePoint< SPACE_DIM > point, bool concreteMove=false)
unsigned AddNode(Node< SPACE_DIM > *pNewNode)
c_vector< double, SPACE_DIM > GetVectorFromAtoB(const c_vector< double, SPACE_DIM > &rLocation1, const c_vector< double, SPACE_DIM > &rLocation2)
virtual void SetUpBoxCollection(double cutOffLength, c_vector< double, 2 *SPACE_DIM > domainSize, int numLocalRows=PETSC_DECIDE, c_vector< bool, SPACE_DIM > isDimPeriodic=zero_vector< bool >(SPACE_DIM))
PeriodicNodesOnlyMesh(c_vector< double, SPACE_DIM > width)
c_vector< double, SPACE_DIM > GetPeriodicWidths() const
double GetWidth(const unsigned &rDimension) const