AbstractTetrahedralElement.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
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 "UblasCustomFunctions.hpp"
00030 
00031 #include "AbstractTetrahedralElement.hpp"
00032 
00033 #include "Exception.hpp"
00034 
00035 #include <cassert>
00036 
00038 // Implementation
00040 
00041 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00042 void AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>::RefreshJacobian(c_matrix<double, SPACE_DIM, ELEMENT_DIM>& rJacobian)
00043 {
00044     if (this->mIsDeleted)
00045     {
00046         EXCEPTION("Attempting to Refresh a deleted element");
00047     }
00048     for (unsigned i=0; i<SPACE_DIM; i++)
00049     {
00050         for (unsigned j=0; j!=ELEMENT_DIM; j++) // Does a j<ELEMENT_DIM without ever having to test j<0U (#186: pointless comparison of unsigned integer with zero)
00051         {
00052             rJacobian(i,j) = this->GetNodeLocation(j+1,i) - this->GetNodeLocation(0,i);
00053         }
00054     }
00055 }
00056 
00057 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00058 AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>::AbstractTetrahedralElement(unsigned index, const std::vector<Node<SPACE_DIM>*>& rNodes)
00059     : AbstractElement<ELEMENT_DIM, SPACE_DIM>(index, rNodes)
00060 {
00061     // Sanity checking
00062     unsigned num_vectices = ELEMENT_DIM+1;
00063 
00064     double det;
00065 
00066     if (SPACE_DIM == ELEMENT_DIM)
00067     {
00068         // This is so we know it's the first time of asking
00069         // Create Jacobian
00070         c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
00071         try
00072         {
00073             CalculateJacobian(jacobian, det);
00074         }
00075         catch (Exception)
00076         {
00077             // if the Jacobian is negative the orientation of the element is probably
00078             // wrong, so swap the last two nodes around.
00079 
00080             this->mNodes[num_vectices-1] = rNodes[num_vectices-2];
00081             this->mNodes[num_vectices-2] = rNodes[num_vectices-1];
00082 
00083             CalculateJacobian(jacobian, det);
00084             // If determinant < 0 then element nodes are listed clockwise.
00085             // We want them anticlockwise.
00086             assert(det > 0.0);
00087         }
00088     }
00089     else 
00090     {
00091         //This is not a full-dimensional element
00092         c_vector<double, SPACE_DIM> weighted_direction;
00093         double det;
00094         CalculateWeightedDirection(weighted_direction, det);
00095     }
00096 
00097 }
00098 
00099 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00100 AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>::AbstractTetrahedralElement(unsigned index)
00101     : AbstractElement<ELEMENT_DIM,SPACE_DIM>(index)
00102 {}
00103 
00104 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00105 void AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>::CalculateJacobian(c_matrix<double, SPACE_DIM, ELEMENT_DIM>& rJacobian, double& rJacobianDeterminant)
00106 {
00107 
00108     assert(ELEMENT_DIM <= SPACE_DIM);
00109     RefreshJacobian(rJacobian);
00110 
00111     {
00112         rJacobianDeterminant = Determinant(rJacobian);
00113         if (rJacobianDeterminant <= DBL_EPSILON)
00114         {
00115             std::stringstream string_stream;
00116             string_stream << "Jacobian determinant is non-positive: "
00117                           << "determinant = " << rJacobianDeterminant
00118                           << " for element " << this->mIndex;
00119             EXCEPTION(string_stream.str());
00120         }
00121     }
00122 }
00123 
00124 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00125 void AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>::CalculateWeightedDirection(c_vector<double, SPACE_DIM>& rWeightedDirection, double& rJacobianDeterminant)
00126 {
00127 
00128     if (ELEMENT_DIM >= SPACE_DIM)
00129     {
00130         assert(ELEMENT_DIM == SPACE_DIM);
00131         EXCEPTION("WeightedDirection undefined for fully dimensional element");
00132     }
00133 
00134     c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
00135     RefreshJacobian(jacobian);
00136 
00137     //At this point we're only dealing with subspace (ELEMENT_DIM < SPACE_DIM) elem
00138     //We assume that the rWeightedDirection vector and rJacobianDeterminant (length of vector)
00139     //are the values from a previous call.
00140 
00141     //This code is only used when ELEMENT_DIM<SPACE_DIM
00142     switch (ELEMENT_DIM)
00143     {
00144         case 0:
00145             //See specialised template for ELEMENT_DIM==0
00146             NEVER_REACHED;
00147             break;
00148         case 1:
00149             // Linear edge in a 2D plane or in 3D
00150 
00151             rWeightedDirection=matrix_column<c_matrix<double,SPACE_DIM,ELEMENT_DIM> >(jacobian, 0);
00152             break;
00153         case 2:
00154             // Surface triangle in a 3d mesh
00155             assert(SPACE_DIM == 3);
00156             rWeightedDirection(0) = -SubDeterminant(jacobian, 0, 2);
00157             rWeightedDirection(1) =  SubDeterminant(jacobian, 1, 2);
00158             rWeightedDirection(2) = -SubDeterminant(jacobian, 2, 2);
00159             break;
00160         default:
00161            ; // Not going to happen
00162     }
00163     rJacobianDeterminant = norm_2(rWeightedDirection);
00164 
00165     if (rJacobianDeterminant < DBL_EPSILON)
00166     {
00167         EXCEPTION("Jacobian determinant is zero");
00168     }
00169 }
00170 
00171 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00172 c_vector<double, SPACE_DIM> AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>::CalculateCentroid() const
00173 {
00174     c_vector<double, SPACE_DIM> centroid = zero_vector<double>(SPACE_DIM);
00175     for (unsigned i=0; i<=ELEMENT_DIM; i++)
00176     {
00177         centroid += this->mNodes[i]->rGetLocation();
00178     }
00179     return centroid/((double)(ELEMENT_DIM + 1));
00180 }
00181 
00182 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00183 void AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>::CalculateInverseJacobian(c_matrix<double, SPACE_DIM, ELEMENT_DIM>& rJacobian, double& rJacobianDeterminant, c_matrix<double, ELEMENT_DIM, SPACE_DIM>& rInverseJacobian)
00184 {
00185     assert(ELEMENT_DIM <= SPACE_DIM);
00186     CalculateJacobian(rJacobian, rJacobianDeterminant);
00187 
00188     // CalculateJacobian should make sure that the determinant is not close to zero (or, in fact, negative)
00189     assert(rJacobianDeterminant > 0.0);
00190     rInverseJacobian = Inverse(rJacobian);
00191 }
00192 
00193 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00194 double AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>::GetVolume(double determinant) const
00195 {
00196     assert(SPACE_DIM == ELEMENT_DIM);
00197 
00198     if (this->mIsDeleted)
00199     {
00200         return 0.0;
00201     }
00202 
00203     double scale_factor = 1.0;
00204 
00205     if (ELEMENT_DIM == 2)
00206     {
00207         scale_factor = 2.0; // both the volume of the canonical triangle is 1/2
00208     }
00209     else if (ELEMENT_DIM == 3)
00210     {
00211         scale_factor = 6.0; // both the volume of the canonical triangle is 1/6
00212     }
00213     return determinant/scale_factor;
00214 }
00215 
00216 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00217 void AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>::GetStiffnessMatrixGlobalIndices(unsigned problemDim, unsigned* pIndices) const
00218 {
00219     for (unsigned local_index=0; local_index<ELEMENT_DIM+1; local_index++)
00220     {
00221         unsigned node = this->GetNodeGlobalIndex(local_index);
00222 
00223         for (unsigned problem_index=0; problem_index<problemDim; problem_index++)
00224         {
00225             //std::cout << local_index*problemDim + problem_index << std::endl;
00226             pIndices[local_index*problemDim + problem_index] = node*problemDim + problem_index;
00227         }
00228     }
00229 }
00230 
00231 
00233 //                  Specialization for 0d elements                  //
00235 
00240 template<unsigned SPACE_DIM>
00241 AbstractTetrahedralElement<0, SPACE_DIM>::AbstractTetrahedralElement(unsigned index, const std::vector<Node<SPACE_DIM>*>& rNodes)
00242     : AbstractElement<0, SPACE_DIM>(index, rNodes)
00243 {
00244     // Sanity checking
00245     //unsigned total_nodes = 1;
00246     assert(this->mNodes.size() == 1);
00247     assert(SPACE_DIM > 0);
00248 
00249     // This is so we know it's the first time of asking
00250     // Create Jacobian
00251     c_vector<double, SPACE_DIM> weighted_direction;
00252     double det;
00253 
00254     CalculateWeightedDirection(weighted_direction, det);
00255 
00256     // If determinant < 0 then element nodes are listed clockwise.
00257     // We want them anticlockwise.
00258     assert(det > 0.0);
00259 }
00260 
00261 template<unsigned SPACE_DIM>
00262 AbstractTetrahedralElement<0, SPACE_DIM>::AbstractTetrahedralElement(unsigned index)
00263     : AbstractElement<0, SPACE_DIM>(index)
00264 {
00265 }
00266 
00267 template<unsigned SPACE_DIM>
00268 void AbstractTetrahedralElement<0, SPACE_DIM>::CalculateWeightedDirection(
00269         c_vector<double, SPACE_DIM>& rWeightedDirection,
00270         double& rJacobianDeterminant)
00271 {
00272     assert(SPACE_DIM > 0);
00273 
00274     // End point of a line
00275     rWeightedDirection = zero_vector<double>(SPACE_DIM);
00276     rWeightedDirection(0) = 1.0;
00277 
00278     rJacobianDeterminant = 1.0;
00279 }
00280 
00281 template<unsigned SPACE_DIM>
00282 c_vector<double, SPACE_DIM> AbstractTetrahedralElement<0, SPACE_DIM>::CalculateCentroid() const
00283 {
00284     c_vector<double, SPACE_DIM> centroid = this->mNodes[0]->rGetLocation();
00285     return centroid;
00286 }
00287 
00288 template<unsigned SPACE_DIM>
00289 void AbstractTetrahedralElement<0, SPACE_DIM>::GetStiffnessMatrixGlobalIndices(unsigned problemDim, unsigned* pIndices) const
00290 {
00291     for (unsigned local_index=0; local_index<1; local_index++)
00292     {
00293         unsigned node = this->GetNodeGlobalIndex(local_index);
00294 
00295         for (unsigned problem_index=0; problem_index<problemDim; problem_index++)
00296         {
00297             //std::cout << local_index*problemDim + problem_index << std::endl;
00298             pIndices[local_index*problemDim + problem_index] = node*problemDim + problem_index;
00299         }
00300     }
00301 }
00302 
00304 // Explicit instantiation
00306 
00307 template class AbstractTetrahedralElement<0,1>;
00308 template class AbstractTetrahedralElement<1,1>;
00309 template class AbstractTetrahedralElement<0,2>;
00310 template class AbstractTetrahedralElement<1,2>;
00311 template class AbstractTetrahedralElement<2,2>;
00312 template class AbstractTetrahedralElement<0,3>;
00313 template class AbstractTetrahedralElement<1,3>;
00314 template class AbstractTetrahedralElement<2,3>;
00315 template class AbstractTetrahedralElement<3,3>;

Generated on Mon Nov 1 12:35:23 2010 for Chaste by  doxygen 1.5.5