AbstractMesh.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include "AbstractMesh.hpp"
00030 
00032 // Implementation
00034 
00035 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00036 AbstractMesh<ELEMENT_DIM, SPACE_DIM>::AbstractMesh()
00037     : mpDistributedVectorFactory(NULL),
00038       mMeshFileBaseName(""),
00039       mMeshChangesDuringSimulation(false)
00040 {
00041 }
00042 
00043 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00044 AbstractMesh<ELEMENT_DIM, SPACE_DIM>::~AbstractMesh()
00045 {
00046     // Iterate over nodes and free the memory
00047     for (unsigned i=0; i<mNodes.size(); i++)
00048     {
00049         delete mNodes[i];
00050     }
00051     if (mpDistributedVectorFactory)
00052     {
00053         delete mpDistributedVectorFactory;
00054     }
00055 }
00056 
00057 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00058 unsigned AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetNumNodes() const
00059 {
00060     return mNodes.size();
00061 }
00062 
00063 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00064 unsigned AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetNumBoundaryNodes() const
00065 {
00066     return mBoundaryNodes.size();
00067 }
00068 
00069 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00070 unsigned AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetNumAllNodes() const
00071 {
00072     return mNodes.size();
00073 }
00074 
00075 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00076 Node<SPACE_DIM>* AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetNode(unsigned index) const
00077 {
00078     unsigned local_index = SolveNodeMapping(index);
00079     return mNodes[local_index];
00080 }
00081 
00082 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00083 Node<SPACE_DIM>* AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetNodeOrHaloNode(unsigned index) const
00084 {
00085     return GetNode(index);
00086 }
00087 
00088 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00089 Node<SPACE_DIM>* AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetNodeFromPrePermutationIndex(unsigned index) const
00090 {
00091     if (mNodesPermutation.empty())
00092     {
00093         return GetNode(index);
00094     }
00095     else
00096     {
00097         return GetNode(mNodesPermutation[index]);
00098     }
00099 }
00100 
00101 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00102 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::ReadNodesPerProcessorFile(const std::string& rNodesPerProcessorFile)
00103 {
00104     NEVER_REACHED;
00105 }
00106 
00107 
00108 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00109 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::SetElementOwnerships()
00110 {
00111     //Does nothing, since an AbstractMesh has no elements
00112 }
00113 
00114 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00115 DistributedVectorFactory* AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetDistributedVectorFactory()
00116 {
00117     if (mpDistributedVectorFactory == NULL)
00118     {
00119         mpDistributedVectorFactory = new DistributedVectorFactory(GetNumNodes());
00120         if (PetscTools::IsParallel())
00121         {
00122             SetElementOwnerships();  //So any parallel implementation with shared mesh has owners set
00123         }
00124     }
00125     return mpDistributedVectorFactory;
00126 }
00127 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00128 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::SetDistributedVectorFactory(DistributedVectorFactory *pFactory)
00129 {
00130     if (mpDistributedVectorFactory)
00131     {
00132         EXCEPTION("Cannot change the mesh's distributed vector factory once it has been set.");
00133     }
00134     if (pFactory->GetNumProcs() != PetscTools::GetNumProcs())
00135     {
00136         EXCEPTION("The distributed vector factory provided to the mesh is for the wrong number of processes.");
00137     }
00138     mpDistributedVectorFactory = pFactory;
00139 }
00140 
00141 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00142 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::PermuteNodes()
00143 {
00144     NEVER_REACHED;
00145 }
00146 
00147 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00148 typename AbstractMesh<ELEMENT_DIM, SPACE_DIM>::BoundaryNodeIterator AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetBoundaryNodeIteratorBegin() const
00149 {
00150     return mBoundaryNodes.begin();
00151 }
00152 
00153 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00154 typename AbstractMesh<ELEMENT_DIM, SPACE_DIM>::BoundaryNodeIterator AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetBoundaryNodeIteratorEnd() const
00155 {
00156     return mBoundaryNodes.end();
00157 }
00158 
00159 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00160 std::string AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetMeshFileBaseName() const
00161 {
00162     if (!IsMeshOnDisk())
00163     {
00164         EXCEPTION("This mesh was not constructed from a file.");
00165     }
00166 
00167     return mMeshFileBaseName;
00168 }
00169 
00170 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00171 bool AbstractMesh<ELEMENT_DIM, SPACE_DIM>::IsMeshOnDisk() const
00172 {
00173     return (mMeshFileBaseName != "");
00174 }
00175 
00176 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00177 const std::vector<unsigned>& AbstractMesh<ELEMENT_DIM, SPACE_DIM>::rGetNodePermutation() const
00178 {
00179     return mNodesPermutation;
00180 }
00181 
00182 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00183 c_vector<double, SPACE_DIM> AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetVectorFromAtoB(
00184     const c_vector<double, SPACE_DIM>& rLocationA, const c_vector<double, SPACE_DIM>& rLocationB)
00185 {
00186     c_vector<double, SPACE_DIM> vector = rLocationB - rLocationA;
00187     return vector;
00188 }
00189 
00190 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00191 double AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetDistanceBetweenNodes(unsigned indexA, unsigned indexB)
00192 {
00193     c_vector<double, SPACE_DIM> vector = GetVectorFromAtoB(mNodes[indexA]->rGetLocation(),
00194                                                            mNodes[indexB]->rGetLocation());
00195     return norm_2(vector);
00196 }
00197 
00198 
00199 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00200 double AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetWidth(const unsigned& rDimension) const
00201 {
00202     assert(rDimension < SPACE_DIM);
00203     return CalculateBoundingBox().GetWidth(rDimension);
00204 }
00205 
00206 
00207 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00208 ChasteCuboid<SPACE_DIM> AbstractMesh<ELEMENT_DIM, SPACE_DIM>::CalculateBoundingBox() const
00209 {
00210     //Set min to DBL_MAX etc.
00211     c_vector<double, SPACE_DIM> minimum_point = scalar_vector<double>(SPACE_DIM, DBL_MAX);
00212     //Set max to -DBL_MAX etc.
00213     c_vector<double, SPACE_DIM> maximum_point=-minimum_point;
00214 
00215     //Iterate through nodes
00217     for (unsigned i=0; i<mNodes.size(); i++)
00218     {
00219         if (!mNodes[i]->IsDeleted())
00220         {
00221             c_vector<double, SPACE_DIM> position  = mNodes[i]->rGetLocation();
00222             //Update max/min
00223             for (unsigned i=0; i<SPACE_DIM; i++)
00224             {
00225                 if (position[i] < minimum_point[i])
00226                 {
00227                     minimum_point[i] = position[i];
00228                 }
00229                 if (position[i] > maximum_point[i])
00230                 {
00231                     maximum_point[i] = position[i];
00232                 }
00233             }
00234         }
00235     }
00236     ChastePoint<SPACE_DIM> min(minimum_point);
00237     ChastePoint<SPACE_DIM> max(maximum_point);
00238 
00239     return ChasteCuboid<SPACE_DIM>(min, max);
00240 }
00241 
00242 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00243 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::Scale(const double xScale, const double yScale, const double zScale)
00244 {
00245     unsigned num_nodes = mNodes.size();
00246 
00247     for (unsigned i=0; i<num_nodes; i++)
00248     {
00249         c_vector<double, SPACE_DIM>& r_location = mNodes[i]->rGetModifiableLocation();
00250         if (SPACE_DIM>=3)
00251         {
00252             r_location[2] *= zScale;
00253         }
00254         if (SPACE_DIM>=2)
00255         {
00256             r_location[1] *= yScale;
00257         }
00258         r_location[0] *= xScale;
00259     }
00260 
00261     RefreshMesh();
00262 }
00263 
00264 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00265 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::Translate(
00266     const double xMovement,
00267     const double yMovement,
00268     const double zMovement)
00269 {
00270     c_vector<double, SPACE_DIM> displacement;
00271 
00272     switch (SPACE_DIM)
00273     {
00274         case 3:
00275             displacement[2] = zMovement;
00276         case 2:
00277             displacement[1] = yMovement;
00278         case 1:
00279             displacement[0] = xMovement;
00280     }
00281 
00282     Translate(displacement);
00283 }
00284 
00285 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00286 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::Translate(const c_vector<double, SPACE_DIM>& rTransVec)
00287 {
00288     unsigned num_nodes = this->GetNumAllNodes();
00289 
00290     for (unsigned i=0; i<num_nodes; i++)
00291     {
00292         c_vector<double, SPACE_DIM>& r_location = this->mNodes[i]->rGetModifiableLocation();
00293         r_location += rTransVec;
00294     }
00295 
00296     RefreshMesh();
00297 }
00298 
00299 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00300 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::Rotate(c_matrix<double , SPACE_DIM, SPACE_DIM> rotationMatrix)
00301 {
00302     unsigned num_nodes = this->GetNumAllNodes();
00303 
00304     for (unsigned i=0; i<num_nodes; i++)
00305     {
00306         c_vector<double, SPACE_DIM>& r_location = this->mNodes[i]->rGetModifiableLocation();
00307         r_location = prod(rotationMatrix, r_location);
00308     }
00309 
00310     RefreshMesh();
00311 }
00312 
00313 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00314 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::Rotate(c_vector<double,3> axis, double angle)
00315 {
00316     assert(SPACE_DIM == 3);
00317     double norm = norm_2(axis);
00318     c_vector<double,3> unit_axis=axis/norm;
00319 
00320     c_matrix<double, SPACE_DIM,SPACE_DIM> rotation_matrix;
00321 
00322     double c = cos(angle);
00323     double s = sin(angle);
00324 
00325     rotation_matrix(0,0) = unit_axis(0)*unit_axis(0)+c*(1-unit_axis(0)*unit_axis(0));
00326     rotation_matrix(0,1) = unit_axis(0)*unit_axis(1)*(1-c) - unit_axis(2)*s;
00327     rotation_matrix(1,0) = unit_axis(0)*unit_axis(1)*(1-c) + unit_axis(2)*s;
00328     rotation_matrix(1,1) = unit_axis(1)*unit_axis(1)+c*(1-unit_axis(1)*unit_axis(1));
00329     rotation_matrix(0,2) = unit_axis(0)*unit_axis(2)*(1-c)+unit_axis(1)*s;
00330     rotation_matrix(1,2) = unit_axis(1)*unit_axis(2)*(1-c)-unit_axis(0)*s;
00331     rotation_matrix(2,0) = unit_axis(0)*unit_axis(2)*(1-c)-unit_axis(1)*s;
00332     rotation_matrix(2,1) = unit_axis(1)*unit_axis(2)*(1-c)+unit_axis(0)*s;
00333     rotation_matrix(2,2) = unit_axis(2)*unit_axis(2)+c*(1-unit_axis(2)*unit_axis(2));
00334 
00335     Rotate(rotation_matrix);
00336 }
00337 
00338 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00339 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::RotateX(const double theta)
00340 {
00341     if (SPACE_DIM != 3)
00342     {
00343         EXCEPTION("This rotation is only valid in 3D");
00344     }
00345     c_matrix<double , SPACE_DIM, SPACE_DIM> x_rotation_matrix=identity_matrix<double>(SPACE_DIM);
00346 
00347     x_rotation_matrix(1,1) = cos(theta);
00348     x_rotation_matrix(1,2) = sin(theta);
00349     x_rotation_matrix(2,1) = -sin(theta);
00350     x_rotation_matrix(2,2) = cos(theta);
00351     Rotate(x_rotation_matrix);
00352 }
00353 
00354 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00355 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::RotateY(const double theta)
00356 {
00357     if (SPACE_DIM != 3)
00358     {
00359         EXCEPTION("This rotation is only valid in 3D");
00360     }
00361     c_matrix<double , SPACE_DIM, SPACE_DIM> y_rotation_matrix=identity_matrix<double>(SPACE_DIM);
00362 
00363     y_rotation_matrix(0,0) = cos(theta);
00364     y_rotation_matrix(0,2) = -sin(theta);
00365     y_rotation_matrix(2,0) = sin(theta);
00366     y_rotation_matrix(2,2) = cos(theta);
00367 
00368 
00369     Rotate(y_rotation_matrix);
00370 }
00371 
00372 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00373 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::RotateZ(const double theta)
00374 {
00375     if (SPACE_DIM < 2)
00376     {
00377         EXCEPTION("This rotation is not valid in less than 2D");
00378     }
00379     c_matrix<double , SPACE_DIM, SPACE_DIM> z_rotation_matrix=identity_matrix<double>(SPACE_DIM);
00380 
00381 
00382     z_rotation_matrix(0,0) = cos(theta);
00383     z_rotation_matrix(0,1) = sin(theta);
00384     z_rotation_matrix(1,0) = -sin(theta);
00385     z_rotation_matrix(1,1) = cos(theta);
00386 
00387     Rotate(z_rotation_matrix);
00388 }
00389 
00390 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00391 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::Rotate(double theta)
00392 {
00393     RotateZ(theta);
00394 }
00395 
00396 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00397 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::RefreshMesh()
00398 {
00399 }
00400 
00401 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00402 bool AbstractMesh<ELEMENT_DIM, SPACE_DIM>::IsMeshChanging() const
00403 {
00404     return mMeshChangesDuringSimulation;
00405 }
00406 
00407 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00408 unsigned AbstractMesh<ELEMENT_DIM, SPACE_DIM>::CalculateMaximumContainingElementsPerProcess() const
00409 {
00410     unsigned max_num=0u;
00411     for (unsigned local_node_index=0; local_node_index<mNodes.size(); local_node_index++)
00412     {
00413         unsigned num=mNodes[local_node_index]->GetNumContainingElements();
00414         if (num>max_num)
00415         {
00416             max_num=num;
00417         }
00418     }
00419     return max_num;
00420 }
00421 
00422 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00423 void AbstractMesh<ELEMENT_DIM, SPACE_DIM>::SetMeshHasChangedSinceLoading()
00424 {
00425     // We just forget what the original file was, which has the desired effect
00426     mMeshFileBaseName = "";
00427 }
00428 
00430 // Explicit instantiation
00432 
00433 template class AbstractMesh<1,1>;
00434 template class AbstractMesh<1,2>;
00435 template class AbstractMesh<1,3>;
00436 template class AbstractMesh<2,2>;
00437 template class AbstractMesh<2,3>;
00438 template class AbstractMesh<3,3>;

Generated on Tue May 31 14:31:48 2011 for Chaste by  doxygen 1.5.5