00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "Element.hpp"
00032
00033
00035
00037
00038
00039 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00040 Element<ELEMENT_DIM, SPACE_DIM>::Element(unsigned index, std::vector<Node<SPACE_DIM>*> nodes)
00041 : AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>(index, nodes)
00042 {
00043 RegisterWithNodes();
00044 }
00045
00046
00053 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00054 Element<ELEMENT_DIM, SPACE_DIM>::Element(const Element &element, const unsigned index)
00055 {
00056 *this = element;
00057 this->mIndex=index;
00058
00059 RegisterWithNodes();
00060 }
00061
00062 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00063 void Element<ELEMENT_DIM, SPACE_DIM>::RegisterWithNodes()
00064 {
00065 for (unsigned i=0; i<this->mNodes.size(); i++)
00066 {
00067 this->mNodes[i]->AddElement(this->mIndex);
00068 }
00069 }
00070
00071 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00072 void Element<ELEMENT_DIM, SPACE_DIM>::MarkAsDeleted()
00073 {
00074 this->mIsDeleted = true;
00075
00076
00077 for (unsigned i=0; i<this->GetNumNodes(); i++)
00078 {
00079 this->mNodes[i]->RemoveElement(this->mIndex);
00080 }
00081 }
00082
00087 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00088 void Element<ELEMENT_DIM, SPACE_DIM>::UpdateNode(const unsigned& rIndex, Node<SPACE_DIM>* pNode)
00089 {
00090 assert(rIndex < this->mNodes.size());
00091
00092
00093 this->mNodes[rIndex]->RemoveElement(this->mIndex);
00094
00095
00096 this->mNodes[rIndex] = pNode;
00097
00098
00099 this->mNodes[rIndex]->AddElement(this->mIndex);
00100 }
00101
00102 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00103 void Element<ELEMENT_DIM, SPACE_DIM>::ResetIndex(unsigned index)
00104 {
00105
00106 for (unsigned i=0; i<this->GetNumNodes(); i++)
00107 {
00108
00109 this->mNodes[i]->RemoveElement(this->mIndex);
00110 }
00111
00112 this->mIndex=index;
00113 RegisterWithNodes();
00114 }
00115
00121 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00122 c_vector<double,SPACE_DIM+1> Element<ELEMENT_DIM, SPACE_DIM>::CalculateCircumsphere()
00123 {
00124
00125
00126
00127
00128
00129
00130
00131
00132 assert(ELEMENT_DIM == SPACE_DIM);
00133 c_vector <double, ELEMENT_DIM> rhs;
00134
00135 c_matrix<double, SPACE_DIM, SPACE_DIM> jacobian;
00136 c_matrix<double, SPACE_DIM, SPACE_DIM> inverse_jacobian;
00137 double jacobian_determinant;
00138
00139 CalculateInverseJacobian(jacobian, jacobian_determinant, inverse_jacobian);
00140
00141 for (unsigned j=0; j<ELEMENT_DIM; j++)
00142 {
00143 double squared_location=0.0;
00144 for (unsigned i=0; i<SPACE_DIM; i++)
00145 {
00146
00147 squared_location += jacobian(i,j)*jacobian(i,j);
00148 }
00149 rhs[j]=squared_location/2.0;
00150 }
00151
00152 c_vector <double, SPACE_DIM> centre;
00153 centre = prod(rhs, inverse_jacobian);
00154 c_vector <double, SPACE_DIM+1> circum;
00155 double squared_radius=0.0;
00156 for (unsigned i=0; i<SPACE_DIM; i++)
00157 {
00158 circum[i] = centre[i] + this->GetNodeLocation(0,i);
00159 squared_radius += centre[i]*centre[i];
00160 }
00161 circum[SPACE_DIM] = squared_radius;
00162
00163 return circum;
00164
00165 }
00166
00167 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00168 c_vector<double,SPACE_DIM+1> Element<ELEMENT_DIM, SPACE_DIM>::CalculateCircumsphere(c_matrix<double, SPACE_DIM, SPACE_DIM>& rJacobian, c_matrix<double, SPACE_DIM, SPACE_DIM>& rInverseJacobian)
00169 {
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179 assert(ELEMENT_DIM == SPACE_DIM);
00180 c_vector <double, ELEMENT_DIM> rhs;
00181
00182 for (unsigned j=0; j<ELEMENT_DIM; j++)
00183 {
00184 double squared_location=0.0;
00185 for (unsigned i=0; i<SPACE_DIM; i++)
00186 {
00187
00188 squared_location += rJacobian(i,j)*rJacobian(i,j);
00189 }
00190 rhs[j]=squared_location/2.0;
00191 }
00192
00193 c_vector <double, SPACE_DIM> centre;
00194 centre = prod(rhs, rInverseJacobian);
00195 c_vector <double, SPACE_DIM+1> circum;
00196 double squared_radius = 0.0;
00197 for (unsigned i=0; i<SPACE_DIM; i++)
00198 {
00199 circum[i] = centre[i] + this->GetNodeLocation(0,i);
00200 squared_radius += centre[i]*centre[i];
00201 }
00202 circum[SPACE_DIM] = squared_radius;
00203
00204 return circum;
00205 }
00206
00207 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00208 double Element<ELEMENT_DIM, SPACE_DIM>::CalculateCircumsphereVolume()
00209 {
00210 c_vector<double, SPACE_DIM+1> circum=CalculateCircumsphere();
00211 if (SPACE_DIM == 1)
00212 {
00213 return 2.0*sqrt(circum[SPACE_DIM]);
00214 }
00215 else if (SPACE_DIM == 2)
00216 {
00217 return M_PI*circum[SPACE_DIM];
00218 }
00219 assert(SPACE_DIM == 3);
00220 return 4.0*M_PI*circum[SPACE_DIM]*sqrt(circum[SPACE_DIM])/3.0;
00221 }
00222
00228 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00229 double Element<ELEMENT_DIM, SPACE_DIM>::CalculateQuality()
00230 {
00231 assert(SPACE_DIM == ELEMENT_DIM);
00232 if (SPACE_DIM == 1)
00233 {
00234 return 1.0;
00235 }
00236
00237 c_matrix<double, SPACE_DIM, SPACE_DIM> jacobian;
00238 double jacobian_determinant;
00239
00240 CalculateJacobian(jacobian, jacobian_determinant);
00241
00242 c_vector<double, SPACE_DIM+1> circum=CalculateCircumsphere();
00243 if (SPACE_DIM == 2)
00244 {
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258 return 2.0*jacobian_determinant/(3.0*sqrt(3)*circum[SPACE_DIM]);
00259 }
00260 assert(SPACE_DIM == 3);
00261
00262
00263
00264
00265
00266
00267
00268
00269 return (3.0*sqrt(3.0)*jacobian_determinant)
00270 /(16.0*circum[SPACE_DIM]*sqrt(circum[SPACE_DIM]));
00271 }
00272
00273 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00274 c_vector<double, SPACE_DIM+1> Element<ELEMENT_DIM, SPACE_DIM>::CalculateInterpolationWeights(ChastePoint<SPACE_DIM> testPoint)
00275 {
00276
00277 assert(ELEMENT_DIM == SPACE_DIM);
00278
00279 c_vector<double, SPACE_DIM+1> weights;
00280
00281 c_vector<double, SPACE_DIM> psi=CalculatePsi(testPoint);
00282
00283
00284 weights[0]=1.0;
00285 for (unsigned i=1; i<=SPACE_DIM; i++)
00286 {
00287 weights[0] -= psi[i-1];
00288 weights[i] = psi[i-1];
00289 }
00290 return weights;
00291 }
00292
00298 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00299 c_vector<double, SPACE_DIM+1> Element<ELEMENT_DIM, SPACE_DIM>::CalculateInterpolationWeightsWithProjection(ChastePoint<SPACE_DIM> testPoint)
00300 {
00301
00302 assert(ELEMENT_DIM == SPACE_DIM);
00303
00304 c_vector<double, SPACE_DIM+1> weights = CalculateInterpolationWeights(testPoint);
00305
00306
00307 bool negative_weight = false;
00308
00309 for(unsigned i=0;i<=SPACE_DIM;i++)
00310 {
00311 if(weights[i] < 0.0)
00312 {
00313 weights[i] = 0.0;
00314
00315 negative_weight = true;
00316 }
00317 }
00318
00319 if(negative_weight == false)
00320 {
00321
00322 return weights;
00323 }
00324
00325
00326
00327
00328 double sum = norm_1 (weights);
00329
00330 assert(sum >= 1.0);
00331
00332 weights = weights/sum;
00333
00334 return weights;
00335
00336 }
00337
00338 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00339 c_vector<double, SPACE_DIM> Element<ELEMENT_DIM, SPACE_DIM>::CalculatePsi(ChastePoint<SPACE_DIM> testPoint)
00340 {
00341
00342 assert(ELEMENT_DIM == SPACE_DIM);
00343
00344
00345 c_vector<double, SPACE_DIM> test_location=testPoint.rGetLocation()-this->GetNodeLocation(0);
00346
00347
00348 c_matrix<double, SPACE_DIM, SPACE_DIM> jacobian;
00349 c_matrix<double, SPACE_DIM, SPACE_DIM> inverse_jacobian;
00350 double jacobian_determinant;
00351
00352 CalculateInverseJacobian(jacobian, jacobian_determinant, inverse_jacobian);
00353
00354 return prod(inverse_jacobian, test_location);
00355 }
00356
00357
00358 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00359 bool Element<ELEMENT_DIM, SPACE_DIM>::IncludesPoint(ChastePoint<SPACE_DIM> testPoint, bool strict)
00360 {
00361
00362 assert(ELEMENT_DIM == SPACE_DIM);
00363
00364 c_vector<double, SPACE_DIM+1> weights=CalculateInterpolationWeights(testPoint);
00365
00366
00367
00368 for (unsigned i=0;i<=SPACE_DIM;i++)
00369 {
00370 if (strict)
00371 {
00372
00373 if (weights[i] <= 2*DBL_EPSILON)
00374 {
00375 return false;
00376 }
00377 }
00378 else
00379 {
00380
00381 if (weights[i] < -2*DBL_EPSILON)
00382 {
00383 return false;
00384 }
00385 }
00386 }
00387 return true;
00388 }
00389
00390 template class Element<1,1>;
00391 template class Element<1,2>;
00392 template class Element<2,2>;
00393 template class Element<2,3>;
00394 template class Element<3,3>;