37#include "AbstractTetrahedralElement.hpp"
46template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
51 EXCEPTION(
"Attempting to Refresh a deleted element");
53 for (
unsigned i=0; i<SPACE_DIM; i++)
55 for (
unsigned j=0; j!=ELEMENT_DIM; j++)
57 rJacobian(i,j) = this->GetNodeLocation(j+1,i) - this->GetNodeLocation(0,i);
62template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
67 unsigned num_vectices = ELEMENT_DIM+1;
71 if (SPACE_DIM == ELEMENT_DIM)
75 c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
85 this->
mNodes[num_vectices-1] = rNodes[num_vectices-2];
86 this->
mNodes[num_vectices-2] = rNodes[num_vectices-1];
97 c_vector<double, SPACE_DIM> weighted_direction;
103template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
108template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
112 assert(ELEMENT_DIM <= SPACE_DIM);
113 RefreshJacobian(rJacobian);
117 if (rJacobianDeterminant <= DBL_EPSILON)
119 std::stringstream message;
120 message <<
"Jacobian determinant is non-positive: "
121 <<
"determinant = " << rJacobianDeterminant
122 <<
" for element " << this->mIndex <<
" (" << ELEMENT_DIM
123 <<
"D element in " << SPACE_DIM <<
"D space). Nodes are at:" << std::endl;
125 for (
unsigned local_node_index=0u; local_node_index != ELEMENT_DIM+1; local_node_index++)
127 c_vector<double, SPACE_DIM> location = this->GetNodeLocation(local_node_index);
128 message <<
"Node " << this->GetNodeGlobalIndex(local_node_index) <<
":\t";
130 for (
unsigned i=0; i<SPACE_DIM; i++)
132 message << location[i];
133 if (i==SPACE_DIM - 1u)
135 message << std::endl;
148template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
152 if (ELEMENT_DIM >= SPACE_DIM)
154 assert(ELEMENT_DIM == SPACE_DIM);
155 EXCEPTION(
"WeightedDirection undefined for fully dimensional element");
158 c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
159 RefreshJacobian(jacobian);
176 rWeightedDirection=matrix_column<c_matrix<double,SPACE_DIM,ELEMENT_DIM> >(jacobian, 0);
180 assert(SPACE_DIM == 3);
188 rJacobianDeterminant = norm_2(rWeightedDirection);
190 if (rJacobianDeterminant < DBL_EPSILON)
192 EXCEPTION(
"Jacobian determinant is zero");
196template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
199 if (ELEMENT_DIM == 1 && SPACE_DIM == 3)
201 EXCEPTION(
"Don't have enough information to calculate a normal vector");
203 c_vector<double, SPACE_DIM> normal;
205 CalculateWeightedDirection(normal, determinant);
206 normal /= determinant;
207 if (ELEMENT_DIM == 1)
210 double x = normal[0];
211 normal[0] = normal[1];
217template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
220 c_vector<double, SPACE_DIM> centroid = zero_vector<double>(SPACE_DIM);
221 for (
unsigned i=0; i<=ELEMENT_DIM; i++)
223 centroid += this->mNodes[i]->rGetLocation();
225 return centroid/((
double)(ELEMENT_DIM + 1));
228template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
231 assert(ELEMENT_DIM <= SPACE_DIM);
232 CalculateJacobian(rJacobian, rJacobianDeterminant);
235 assert(rJacobianDeterminant > 0.0);
236 rInverseJacobian =
Inverse(rJacobian);
239template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
243 if (this->mIsDeleted)
248 double scale_factor = 1.0;
250 if (ELEMENT_DIM == 2)
254 else if (ELEMENT_DIM == 3)
258 return determinant/scale_factor;
261template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
264 for (
unsigned local_index=0; local_index<ELEMENT_DIM+1; local_index++)
266 unsigned node = this->GetNodeGlobalIndex(local_index);
268 for (
unsigned problem_index=0; problem_index<problemDim; problem_index++)
271 pIndices[local_index*problemDim + problem_index] = node*problemDim + problem_index;
285template<
unsigned SPACE_DIM>
291 assert(this->
mNodes.size() == 1);
292 assert(SPACE_DIM > 0);
296 c_vector<double, SPACE_DIM> weighted_direction;
306template<
unsigned SPACE_DIM>
312template<
unsigned SPACE_DIM>
314 c_vector<double, SPACE_DIM>& rWeightedDirection,
315 double& rJacobianDeterminant)
317 assert(SPACE_DIM > 0);
320 rWeightedDirection = zero_vector<double>(SPACE_DIM);
321 rWeightedDirection(0) = 1.0;
323 rJacobianDeterminant = 1.0;
326template<
unsigned SPACE_DIM>
329 assert(SPACE_DIM > 0);
332 c_vector<double, SPACE_DIM> normal = zero_vector<double>(SPACE_DIM);
337template<
unsigned SPACE_DIM>
340 c_vector<double, SPACE_DIM> centroid = this->mNodes[0]->rGetLocation();
344template<
unsigned SPACE_DIM>
347 for (
unsigned local_index=0; local_index<1; local_index++)
349 unsigned node = this->GetNodeGlobalIndex(local_index);
351 for (
unsigned problem_index=0; problem_index<problemDim; problem_index++)
354 pIndices[local_index*problemDim + problem_index] = node*problemDim + problem_index;
#define EXCEPTION(message)
T Determinant(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
T SubDeterminant(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM, const unsigned missrow, const unsigned misscol)
boost::numeric::ublas::c_matrix< T, 1, 1 > Inverse(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
std::vector< Node< SPACE_DIM > * > mNodes
void RefreshJacobian(c_matrix< double, SPACE_DIM, ELEMENT_DIM > &rJacobian)
void CalculateInverseJacobian(c_matrix< double, SPACE_DIM, ELEMENT_DIM > &rJacobian, double &rJacobianDeterminant, c_matrix< double, ELEMENT_DIM, SPACE_DIM > &rInverseJacobian)
double GetVolume(double determinant) const
void CalculateWeightedDirection(c_vector< double, SPACE_DIM > &rWeightedDirection, double &rJacobianDeterminant)
AbstractTetrahedralElement(unsigned index, const std::vector< Node< SPACE_DIM > * > &rNodes)
c_vector< double, SPACE_DIM > CalculateCentroid() const
void GetStiffnessMatrixGlobalIndices(unsigned problemDim, unsigned *pIndices) const
void CalculateJacobian(c_matrix< double, SPACE_DIM, ELEMENT_DIM > &rJacobian, double &rJacobianDeterminant)
c_vector< double, SPACE_DIM > CalculateNormal()