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 #ifndef _ABSTRACTTETRAHEDRALELEMENT_HPP_
00031 #define _ABSTRACTTETRAHEDRALELEMENT_HPP_
00032
00033 #include "AbstractElement.hpp"
00034
00038 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00039 class AbstractTetrahedralElement : public AbstractElement<ELEMENT_DIM,SPACE_DIM>
00040 {
00041 protected:
00042
00043 void RefreshJacobian(c_matrix<double, SPACE_DIM, SPACE_DIM>& rJacobian)
00044 {
00045 if (this->mIsDeleted)
00046 {
00047 EXCEPTION("Attempting to Refresh a deleted element");
00048 }
00049 for (unsigned i=0; i<SPACE_DIM; i++)
00050 {
00051 for (unsigned j=0; j!=ELEMENT_DIM; j++)
00052 {
00053 rJacobian(i,j) = this->GetNodeLocation(j+1,i) - this->GetNodeLocation(0,i);
00054 }
00055 }
00056 }
00057
00058 public:
00060 AbstractTetrahedralElement(unsigned index, const std::vector<Node<SPACE_DIM>*>& rNodes);
00061
00066 AbstractTetrahedralElement(unsigned index=INDEX_IS_NOT_USED)
00067 : AbstractElement<ELEMENT_DIM,SPACE_DIM>(index)
00068 {}
00069
00070 virtual ~AbstractTetrahedralElement()
00071 {}
00072
00073 void ZeroJacobianDeterminant(void);
00074 void ZeroWeightedDirection(void);
00075
00076
00077 c_vector<double, SPACE_DIM> CalculateCentroid() const
00078 {
00079 c_vector<double, SPACE_DIM> centroid=zero_vector<double>(SPACE_DIM);
00080 for (unsigned i=0; i<=ELEMENT_DIM; i++)
00081 {
00082 centroid += this->mNodes[i]->rGetLocation();
00083 }
00084 return centroid/((double)(ELEMENT_DIM + 1));
00085 }
00086
00088 void CalculateJacobian(c_matrix<double, SPACE_DIM, SPACE_DIM>& rJacobian, double &rJacobianDeterminant, bool concreteMove=true);
00089 void CalculateWeightedDirection(c_vector<double, SPACE_DIM>& rWeightedDirection, double &rJacobianDeterminant, bool concreteMove=true);
00090
00091
00092 void CalculateInverseJacobian(c_matrix<double, SPACE_DIM, SPACE_DIM>& rJacobian, double &rJacobianDeterminant, c_matrix<double, SPACE_DIM, SPACE_DIM>& rInverseJacobian)
00093 {
00094 assert(ELEMENT_DIM==SPACE_DIM);
00095 CalculateJacobian(rJacobian, rJacobianDeterminant);
00096
00097 assert(rJacobianDeterminant > 0.0);
00098 rInverseJacobian = Inverse(rJacobian);
00099 }
00100
00101
00103
00104 double GetVolume(void)
00105 {
00106 assert(SPACE_DIM == ELEMENT_DIM);
00107
00108 if (this->mIsDeleted)
00109 {
00110 return 0.0;
00111 }
00112
00113
00115 c_matrix<double, SPACE_DIM, SPACE_DIM> jacobian;
00116 double determinant;
00117
00118 CalculateJacobian(jacobian, determinant);
00119 double scale_factor = 1.0;
00120
00121 if (ELEMENT_DIM == 2)
00122 {
00123 scale_factor = 2.0;
00124 }
00125 else if (ELEMENT_DIM == 3)
00126 {
00127 scale_factor= 6.0;
00128 }
00129 return determinant/scale_factor;
00130 }
00131
00140 void GetStiffnessMatrixGlobalIndices(unsigned problemDim, unsigned* pIndices) const
00141 {
00142 for (unsigned local_index=0; local_index<ELEMENT_DIM+1; local_index++)
00143 {
00144 unsigned node = this->GetNodeGlobalIndex(local_index);
00145
00146 for (unsigned problem_index=0; problem_index<problemDim; problem_index++)
00147 {
00148
00149 pIndices[local_index*problemDim + problem_index] = node*problemDim + problem_index;
00150 }
00151 }
00152 }
00153 };
00154
00155
00156
00157 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00158 AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>::AbstractTetrahedralElement(unsigned index,
00159 const std::vector<Node<SPACE_DIM>*>& rNodes)
00160 : AbstractElement<ELEMENT_DIM, SPACE_DIM>(index, rNodes)
00161 {
00162
00163 unsigned total_nodes = ELEMENT_DIM+1;
00164 assert(this->mNodes.size() == total_nodes);
00165
00166
00167
00169 c_matrix<double, SPACE_DIM, SPACE_DIM> jacobian;
00170 c_vector<double, SPACE_DIM> weighted_direction;
00171 double det;
00172
00173 if (SPACE_DIM == ELEMENT_DIM)
00174 {
00175 try
00176 {
00177 CalculateJacobian(jacobian, det);
00178 }
00179 catch (Exception)
00180 {
00181
00182
00183
00184 this->mNodes[total_nodes-1] = rNodes[total_nodes-2];
00185 this->mNodes[total_nodes-2] = rNodes[total_nodes-1];
00186
00187 CalculateJacobian(jacobian, det);
00188 }
00189 }
00190 else
00191 {
00192 CalculateWeightedDirection(weighted_direction, det);
00193 }
00194
00195
00196
00197 assert(det > 0.0);
00198 }
00199
00200
00201
00202
00203
00204 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00205 void AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>::CalculateJacobian(c_matrix<double, SPACE_DIM, SPACE_DIM>& rJacobian, double &rJacobianDeterminant, bool concreteMove)
00206 {
00207
00208 assert(ELEMENT_DIM == SPACE_DIM);
00209 RefreshJacobian(rJacobian);
00210
00211 {
00212 rJacobianDeterminant = Determinant(rJacobian);
00213 if (rJacobianDeterminant <= DBL_EPSILON)
00214 {
00215 std::stringstream string_stream;
00216 string_stream << "Jacobian determinant is non-positive: "
00217 << "determinant = " << rJacobianDeterminant
00218 << " for element " << this->mIndex;
00219 EXCEPTION(string_stream.str());
00220 }
00221 }
00222 }
00223
00224 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00225 void AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>::CalculateWeightedDirection(c_vector<double, SPACE_DIM>& rWeightedDirection, double &rJacobianDeterminant, bool concreteMove)
00226 {
00227 if(ELEMENT_DIM >= SPACE_DIM)
00228 {
00229 assert(ELEMENT_DIM == SPACE_DIM);
00230 EXCEPTION("WeightedDirection undefined for fully dimensional element");
00231 }
00232
00233 c_matrix<double, SPACE_DIM, SPACE_DIM> jacobian;
00234 RefreshJacobian(jacobian);
00235
00236
00237
00238
00239
00240 c_vector<double, SPACE_DIM> weighted_direction;
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250 switch (ELEMENT_DIM)
00251 {
00252 case 0:
00253
00254 weighted_direction(0)=1.0;
00255 if (SPACE_DIM == 2)
00256 {
00257 weighted_direction(1)=0.0;
00258 }
00259 break;
00260 case 1:
00261
00262
00263 weighted_direction=matrix_column<c_matrix<double,SPACE_DIM,SPACE_DIM> >(jacobian,0);
00264 break;
00265 case 2:
00266
00267 assert(SPACE_DIM == 3);
00268 weighted_direction(0)=-SubDeterminant(jacobian,0,2);
00269 weighted_direction(1)= SubDeterminant(jacobian,1,2);
00270 weighted_direction(2)=-SubDeterminant(jacobian,2,2);
00271 break;
00272 default:
00273 ;
00274 }
00275 double jacobian_determinant = norm_2(weighted_direction);
00276 if (jacobian_determinant < DBL_EPSILON)
00277 {
00278 EXCEPTION("Jacobian determinant is zero");
00279 }
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289 rJacobianDeterminant = jacobian_determinant;
00290 rWeightedDirection = weighted_direction;
00291
00292 }
00293
00294 #endif //_ABSTRACTTETRAHEDRALELEMENT_HPP_