Element.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 
00030 
00031 #include "Element.hpp"
00032 
00033 #include <cassert>
00034 
00035 
00037 // Implementation
00039 
00040 
00041 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00042 Element<ELEMENT_DIM, SPACE_DIM>::Element(unsigned index, const std::vector<Node<SPACE_DIM>*>& rNodes)
00043     : AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>(index, rNodes)
00044 {
00045     RegisterWithNodes();
00046 }
00047 
00048 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00049 Element<ELEMENT_DIM, SPACE_DIM>::Element(const Element& rElement, const unsigned index)
00050 {
00051     *this = rElement;
00052     this->mIndex = index;
00053 
00054     RegisterWithNodes();
00055 }
00056 
00057 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00058 void Element<ELEMENT_DIM, SPACE_DIM>::RegisterWithNodes()
00059 {
00060     for (unsigned i=0; i<this->mNodes.size(); i++)
00061     {
00062         this->mNodes[i]->AddElement(this->mIndex);
00063     }
00064 }
00065 
00066 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00067 void Element<ELEMENT_DIM, SPACE_DIM>::MarkAsDeleted()
00068 {
00069     this->mIsDeleted = true;
00070     // Update nodes in this element so they know they are not contained by us
00071     for (unsigned i=0; i<this->GetNumNodes(); i++)
00072     {
00073         this->mNodes[i]->RemoveElement(this->mIndex);
00074     }
00075 }
00076 
00081 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00082 void Element<ELEMENT_DIM, SPACE_DIM>::UpdateNode(const unsigned& rIndex, Node<SPACE_DIM>* pNode)
00083 {
00084     assert(rIndex < this->mNodes.size());
00085 
00086     // Remove it from the node at this location
00087     this->mNodes[rIndex]->RemoveElement(this->mIndex);
00088 
00089     // Update the node at this location
00090     this->mNodes[rIndex] = pNode;
00091 
00092     // Add element to this node
00093     this->mNodes[rIndex]->AddElement(this->mIndex);
00094 }
00095 
00096 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00097 void Element<ELEMENT_DIM, SPACE_DIM>::ResetIndex(unsigned index)
00098 {
00099     //std::cout << "ResetIndex - removing nodes.\n" << std::flush;
00100     for (unsigned i=0; i<this->GetNumNodes(); i++)
00101     {
00102        //std::cout << "Node " << this->mNodes[i]->GetIndex() << " element "<< this->mIndex << std::flush;
00103        this->mNodes[i]->RemoveElement(this->mIndex);
00104     }
00105     //std::cout << "\nResetIndex - done.\n" << std::flush;
00106     this->mIndex=index;
00107     RegisterWithNodes();
00108 }
00109 
00110 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00111 c_vector<double,SPACE_DIM+1> Element<ELEMENT_DIM, SPACE_DIM>::CalculateCircumsphere(c_matrix<double, SPACE_DIM, ELEMENT_DIM>& rJacobian, c_matrix<double, ELEMENT_DIM, SPACE_DIM>& rInverseJacobian)
00112 {
00113     /*Assuming that x0,y0.. is at the origin then we need to solve
00114      *
00115      * ( 2x1 2y1 2z1  ) (x)    (x1^2+y1^2+z1^2)
00116      * ( 2x2 2y2 2z2  ) (y)    (x2^2+y2^2+z2^2)
00117      * ( 2x3 2y3 2z3  ) (z)    (x3^2+y3^2+z3^2)
00118      * where (x,y,z) is the circumcentre
00119      *
00120      */
00121 
00122     assert(ELEMENT_DIM == SPACE_DIM);
00123     c_vector<double, ELEMENT_DIM> rhs;
00124 
00125     for (unsigned j=0; j<ELEMENT_DIM; j++)
00126     {
00127         double squared_location=0.0;
00128         for (unsigned i=0; i<SPACE_DIM; i++)
00129         {
00130             //mJacobian(i,j) is the i-th component of j-th vertex (relative to vertex 0)
00131             squared_location += rJacobian(i,j)*rJacobian(i,j);
00132         }
00133         rhs[j]=squared_location/2.0;
00134     }
00135 
00136     c_vector<double, SPACE_DIM> centre;
00137     centre = prod(rhs, rInverseJacobian);
00138     c_vector<double, SPACE_DIM+1> circum;
00139     double squared_radius = 0.0;
00140     for (unsigned i=0; i<SPACE_DIM; i++)
00141     {
00142         circum[i] = centre[i] + this->GetNodeLocation(0,i);
00143         squared_radius += centre[i]*centre[i];
00144     }
00145     circum[SPACE_DIM] = squared_radius;
00146 
00147     return circum;
00148 }
00149 
00155 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00156 double Element<ELEMENT_DIM, SPACE_DIM>::CalculateQuality()
00157 {
00158     assert(SPACE_DIM == ELEMENT_DIM);
00159     if (SPACE_DIM == 1)
00160     {
00161         return 1.0;
00162     }
00163 
00164     c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
00165     c_matrix<double, ELEMENT_DIM, SPACE_DIM> jacobian_inverse;
00166     double jacobian_determinant;
00167 
00168     this->CalculateInverseJacobian(jacobian, jacobian_determinant, jacobian_inverse);
00169 
00170     c_vector<double, SPACE_DIM+1> circum=CalculateCircumsphere(jacobian, jacobian_inverse);
00171     if (SPACE_DIM == 2)
00172     {
00173         /* Want Q=(Area_Tri / Area_Cir) / (Area_Equilateral_Tri / Area_Equilateral_Cir)
00174          * Area_Tri = |Jacobian| /2
00175          * Area_Cir = Pi * r^2
00176          * Area_Eq_Tri = (3*sqrt(3)/4)*R^2
00177          * Area_Eq_Tri = Pi * R^2
00178          * Q= (2*|Jacobian|)/3*sqrt(3)*r^2)
00179          */
00180         return 2.0*jacobian_determinant/(3.0*sqrt(3.0)*circum[SPACE_DIM]);
00181     }
00182     assert(SPACE_DIM == 3);
00183     /* Want Q=(Vol_Tet / Vol_CirS) / (Vol_Plat_Tet / Vol_Plat_CirS)
00184       *  Vol_Tet  = |Jacobian| /6
00185       *  Vol_CirS = 4*Pi*r^3/3
00186       *  Vol_Plat_Tet  = 8*sqrt(3)*R^3/27
00187       *  Vol_Plat_CirS = 4*Pi*R^3/3
00188      * Q= 3*sqrt(3)*|Jacobian|/ (16*r^3)
00189       */
00190 
00191     return (3.0*sqrt(3.0)*jacobian_determinant)
00192            /(16.0*circum[SPACE_DIM]*sqrt(circum[SPACE_DIM]));
00193 }
00194 
00195 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00196 c_vector<double, SPACE_DIM+1> Element<ELEMENT_DIM, SPACE_DIM>::CalculateInterpolationWeights(const ChastePoint<SPACE_DIM>& rTestPoint)
00197 {
00198     //Can only test if it's a tetrahedal mesh in 3d, triangles in 2d...
00199     assert(ELEMENT_DIM == SPACE_DIM);
00200 
00201     c_vector<double, SPACE_DIM+1> weights;
00202 
00203     c_vector<double, SPACE_DIM> xi=CalculateXi(rTestPoint);
00204 
00205     //Copy 3 weights and compute the fourth weight
00206     weights[0]=1.0;
00207     for (unsigned i=1; i<=SPACE_DIM; i++)
00208     {
00209         weights[0] -= xi[i-1];
00210         weights[i] = xi[i-1];
00211     }
00212     return weights;
00213 }
00214 
00220 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00221 c_vector<double, SPACE_DIM+1> Element<ELEMENT_DIM, SPACE_DIM>::CalculateInterpolationWeightsWithProjection(const ChastePoint<SPACE_DIM>& rTestPoint)
00222 {
00223     //Can only test if it's a tetrahedal mesh in 3d, triangles in 2d...
00224     assert(ELEMENT_DIM == SPACE_DIM);
00225 
00226     c_vector<double, SPACE_DIM+1> weights = CalculateInterpolationWeights(rTestPoint);
00227 
00228     // Check for negative weights and set them to zero.
00229     bool negative_weight = false;
00230 
00231     for (unsigned i=0; i<=SPACE_DIM; i++)
00232     {
00233         if (weights[i] < 0.0)
00234         {
00235             weights[i] = 0.0;
00236 
00237             negative_weight = true;
00238         }
00239     }
00240 
00241     if (negative_weight == false)
00242     {
00243         // If there are no negative weights, there is nothing to do.
00244         return weights;
00245     }
00246 
00247     // Renormalise so that all weights add to 1.0.
00248 
00249     // Note that all elements of weights are now non-negative and so the l1-norm (sum of magnitudes) is equivalent to the sum of the elements of the vector
00250     double sum = norm_1 (weights);
00251 
00252     //l1-norm ought to be above 1 (because we scrubbed negative weights)
00253     //However, if we scrubbed weights that were the size of the machine precision then we might be close to one (even less than 1).
00254     assert( sum + DBL_EPSILON >= 1.0);
00255 
00256     //We might skip this division when sum ~= 1
00257     weights = weights/sum;
00258 
00259     return weights;
00260 }
00261 
00262 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00263 c_vector<double, SPACE_DIM> Element<ELEMENT_DIM, SPACE_DIM>::CalculateXi(const ChastePoint<SPACE_DIM>& rTestPoint)
00264 {
00265     //Can only test if it's a tetrahedal mesh in 3d, triangles in 2d...
00266     assert(ELEMENT_DIM == SPACE_DIM);
00267 
00268     //Find the location with respect to node 0
00269 // todo: #1361 ComputeContainingElements and related methods, and methods called by that down to 
00270 // here, should really take in const c_vector& rather than ChastePoints.
00271     ChastePoint<SPACE_DIM> copy = rTestPoint; // as rGetLocation in line below is not a const method
00272     c_vector<double, SPACE_DIM> test_location=copy.rGetLocation()-this->GetNodeLocation(0);
00273 
00274     //Multiply by inverse Jacobian
00275     c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
00276     c_matrix<double, ELEMENT_DIM, SPACE_DIM> inverse_jacobian;
00277     double jacobian_determinant;
00278 
00280     this->CalculateInverseJacobian(jacobian, jacobian_determinant, inverse_jacobian);
00281 
00282     return prod(inverse_jacobian, test_location);
00283 }
00284 
00285 
00286 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00287 bool Element<ELEMENT_DIM, SPACE_DIM>::IncludesPoint(const ChastePoint<SPACE_DIM>& rTestPoint, bool strict)
00288 {
00289     //Can only test if it's a tetrahedal mesh in 3d, triangles in 2d...
00290     assert(ELEMENT_DIM == SPACE_DIM);
00291 
00292     c_vector<double, SPACE_DIM+1> weights=CalculateInterpolationWeights(rTestPoint);
00293 
00294     //If the point is in the simplex then all the weights should be positive
00295 
00296     for (unsigned i=0; i<=SPACE_DIM; i++)
00297     {
00298         if (strict)
00299         {
00300             //Points can't be close to a face
00301             if (weights[i] <= 2*DBL_EPSILON)
00302             {
00303                 return false;
00304             }
00305         }
00306         else
00307         {
00308             //Allow point to be close to a face
00309             if (weights[i] < -2*DBL_EPSILON)
00310             {
00311                 return false;
00312             }
00313         }
00314     }
00315     return true;
00316 }
00317 
00319 // Explicit instantiation
00321 
00322 template class Element<1,1>;
00323 template class Element<1,2>;
00324 template class Element<1,3>;
00325 template class Element<2,2>;
00326 template class Element<2,3>;
00327 template class Element<3,3>;

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