Chaste  Release::3.4
Element.cpp
1 /*
2 
3 Copyright (c) 2005-2016, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "Element.hpp"
37 
38 #include <cfloat>
39 #include <cassert>
40 
42 // Implementation
44 
45 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
46 Element<ELEMENT_DIM, SPACE_DIM>::Element(unsigned index, const std::vector<Node<SPACE_DIM>*>& rNodes, bool registerWithNodes)
47  : AbstractTetrahedralElement<ELEMENT_DIM, SPACE_DIM>(index, rNodes)
48 {
49  if (registerWithNodes)
50  {
52  }
53 }
54 
55 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
56 Element<ELEMENT_DIM, SPACE_DIM>::Element(const Element& rElement, const unsigned index)
57 {
58  *this = rElement;
59  this->mIndex = index;
60 
61  RegisterWithNodes();
62 }
63 
64 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
66 {
67  for (unsigned i=0; i<this->mNodes.size(); i++)
68  {
69  this->mNodes[i]->AddElement(this->mIndex);
70  }
71 }
72 
73 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
75 {
76  this->mIsDeleted = true;
77  // Update nodes in this element so they know they are not contained by us
78  for (unsigned i=0; i<this->GetNumNodes(); i++)
79  {
80  this->mNodes[i]->RemoveElement(this->mIndex);
81  }
82 }
83 
88 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
90 {
91  assert(rIndex < this->mNodes.size());
92 
93  // Remove it from the node at this location
94  this->mNodes[rIndex]->RemoveElement(this->mIndex);
95 
96  // Update the node at this location
97  this->mNodes[rIndex] = pNode;
98 
99  // Add element to this node
100  this->mNodes[rIndex]->AddElement(this->mIndex);
101 }
102 
103 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
105 {
106  //std::cout << "ResetIndex - removing nodes.\n" << std::flush;
107  for (unsigned i=0; i<this->GetNumNodes(); i++)
108  {
109  //std::cout << "Node " << this->mNodes[i]->GetIndex() << " element "<< this->mIndex << std::flush;
110  this->mNodes[i]->RemoveElement(this->mIndex);
111  }
112  //std::cout << "\nResetIndex - done.\n" << std::flush;
113  this->mIndex=index;
114  RegisterWithNodes();
115 }
116 
117 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
118 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)
119 {
120  /*Assuming that x0,y0.. is at the origin then we need to solve
121  *
122  * ( 2x1 2y1 2z1 ) (x) (x1^2+y1^2+z1^2)
123  * ( 2x2 2y2 2z2 ) (y) (x2^2+y2^2+z2^2)
124  * ( 2x3 2y3 2z3 ) (z) (x3^2+y3^2+z3^2)
125  * where (x,y,z) is the circumcentre
126  *
127  */
128 
129  assert(ELEMENT_DIM == SPACE_DIM);
130  c_vector<double, ELEMENT_DIM> rhs;
131 
132  for (unsigned j=0; j<ELEMENT_DIM; j++)
133  {
134  double squared_location=0.0;
135  for (unsigned i=0; i<SPACE_DIM; i++)
136  {
137  //mJacobian(i,j) is the i-th component of j-th vertex (relative to vertex 0)
138  squared_location += rJacobian(i,j)*rJacobian(i,j);
139  }
140  rhs[j]=squared_location/2.0;
141  }
142 
143  c_vector<double, SPACE_DIM> centre;
144  centre = prod(rhs, rInverseJacobian);
145  c_vector<double, SPACE_DIM+1> circum;
146  double squared_radius = 0.0;
147  for (unsigned i=0; i<SPACE_DIM; i++)
148  {
149  circum[i] = centre[i] + this->GetNodeLocation(0,i);
150  squared_radius += centre[i]*centre[i];
151  }
152  circum[SPACE_DIM] = squared_radius;
153 
154  return circum;
155 }
156 
162 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
164 {
165  assert(SPACE_DIM == ELEMENT_DIM);
166  if (SPACE_DIM == 1)
167  {
168  return 1.0;
169  }
170 
171  c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
172  c_matrix<double, ELEMENT_DIM, SPACE_DIM> jacobian_inverse;
173  double jacobian_determinant;
174 
175  this->CalculateInverseJacobian(jacobian, jacobian_determinant, jacobian_inverse);
176 
177  c_vector<double, SPACE_DIM+1> circum=CalculateCircumsphere(jacobian, jacobian_inverse);
178  if (SPACE_DIM == 2)
179  {
180  /* Want Q=(Area_Tri / Area_Cir) / (Area_Equilateral_Tri / Area_Equilateral_Cir)
181  * Area_Tri = |Jacobian| /2
182  * Area_Cir = Pi * r^2
183  * Area_Eq_Tri = (3*sqrt(3.0)/4)*R^2
184  * Area_Eq_Tri = Pi * R^2
185  * Q= (2*|Jacobian|)/3*sqrt(3.0)*r^2)
186  */
187  return 2.0*jacobian_determinant/(3.0*sqrt(3.0)*circum[SPACE_DIM]);
188  }
189  assert(SPACE_DIM == 3);
190  /* Want Q=(Vol_Tet / Vol_CirS) / (Vol_Plat_Tet / Vol_Plat_CirS)
191  * Vol_Tet = |Jacobian| /6
192  * Vol_CirS = 4*Pi*r^3/3
193  * Vol_Plat_Tet = 8*sqrt(3.0)*R^3/27
194  * Vol_Plat_CirS = 4*Pi*R^3/3
195  * Q= 3*sqrt(3.0)*|Jacobian|/ (16*r^3)
196  */
197 
198  return (3.0*sqrt(3.0)*jacobian_determinant)
199  /(16.0*circum[SPACE_DIM]*sqrt(circum[SPACE_DIM]));
200 }
201 
202 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
204 {
205  c_vector <double, 2> min_max;
206  min_max[0] = DBL_MAX; //Min initialised to very large
207  min_max[1] = 0.0; //Max initialised to zero
208  for (unsigned i=0; i<=ELEMENT_DIM; i++)
209  {
210  c_vector<double, SPACE_DIM> loc_i = this->GetNodeLocation(i);
211  for (unsigned j=i+1; j<=ELEMENT_DIM; j++)
212  {
213  double length = norm_2(this->GetNodeLocation(j) - loc_i);
214  if (length < min_max[0])
215  {
216  min_max[0] = length;
217  }
218  if (length > min_max[1])
219  {
220  min_max[1] = length;
221  }
222  }
223  }
224  return min_max;
225 }
226 
227 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
229 {
230  // Can only test if it's a tetrahedral mesh in 3d, triangles in 2d...
231  assert(ELEMENT_DIM == SPACE_DIM);
232 
233  c_vector<double, SPACE_DIM+1> weights;
234 
235  c_vector<double, SPACE_DIM> xi=CalculateXi(rTestPoint);
236 
237  // Copy 3 weights and compute the fourth weight
238  weights[0]=1.0;
239  for (unsigned i=1; i<=SPACE_DIM; i++)
240  {
241  weights[0] -= xi[i-1];
242  weights[i] = xi[i-1];
243  }
244  return weights;
245 }
246 
247 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
249 {
250  //Can only test if it's a tetrahedral mesh in 3d, triangles in 2d...
251  assert(ELEMENT_DIM == SPACE_DIM);
252 
253  c_vector<double, SPACE_DIM+1> weights = CalculateInterpolationWeights(rTestPoint);
254 
255  // Check for negative weights and set them to zero.
256  bool negative_weight = false;
257 
258  for (unsigned i=0; i<=SPACE_DIM; i++)
259  {
260  if (weights[i] < 0.0)
261  {
262  weights[i] = 0.0;
263 
264  negative_weight = true;
265  }
266  }
267 
268  if (negative_weight == false)
269  {
270  // If there are no negative weights, there is nothing to do.
271  return weights;
272  }
273 
274  // Renormalise so that all weights add to 1.0.
275 
276  // 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
277  double sum = norm_1 (weights);
278 
279  //l1-norm ought to be above 1 (because we scrubbed negative weights)
280  //However, if we scrubbed weights that were the size of the machine precision then we might be close to one (even less than 1).
281  assert( sum + DBL_EPSILON >= 1.0);
282 
283  //We might skip this division when sum ~= 1
284  weights = weights/sum;
285 
286  return weights;
287 }
288 
289 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
290 c_vector<double, SPACE_DIM> Element<ELEMENT_DIM, SPACE_DIM>::CalculateXi(const ChastePoint<SPACE_DIM>& rTestPoint)
291 {
292  //Can only test if it's a tetrahedral mesh in 3d, triangles in 2d...
293  assert(ELEMENT_DIM == SPACE_DIM);
294 
295  // Find the location with respect to node 0
298  c_vector<double, SPACE_DIM> test_location=rTestPoint.rGetLocation()-this->GetNodeLocation(0);
299 
300  //Multiply by inverse Jacobian
301  c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
302  c_matrix<double, ELEMENT_DIM, SPACE_DIM> inverse_jacobian;
303  double jacobian_determinant;
304 
306  this->CalculateInverseJacobian(jacobian, jacobian_determinant, inverse_jacobian);
307 
308  return prod(inverse_jacobian, test_location);
309 }
310 
311 
312 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
314 {
315  // Can only test if it's a tetrahedral mesh in 3d, triangles in 2d...
316  assert(ELEMENT_DIM == SPACE_DIM);
317 
318  c_vector<double, SPACE_DIM+1> weights=CalculateInterpolationWeights(rTestPoint);
319 
320  // If the point is in the simplex then all the weights should be positive.
321 
322  for (unsigned i=0; i<=SPACE_DIM; i++)
323  {
324  if (strict)
325  {
326  // Points can't be close to a face
327  if (weights[i] <= 2*DBL_EPSILON)
328  {
329  return false;
330  }
331  }
332  else
333  {
334  // Allow point to be close to a face
335  if (weights[i] < -2*DBL_EPSILON)
336  {
337  return false;
338  }
339  }
340  }
341  return true;
342 }
343 
345 // Explicit instantiation
347 
348 template class Element<1,1>;
349 template class Element<1,2>;
350 template class Element<1,3>;
351 template class Element<2,2>;
352 template class Element<2,3>;
353 template class Element<3,3>;
c_vector< double, SPACE_DIM+1 > CalculateCircumsphere(c_matrix< double, SPACE_DIM, ELEMENT_DIM > &rJacobian, c_matrix< double, ELEMENT_DIM, SPACE_DIM > &rInverseJacobian)
Definition: Element.cpp:118
Definition: Node.hpp:58
c_vector< double, DIM > & rGetLocation()
Definition: ChastePoint.cpp:76
double CalculateQuality()
Definition: Element.cpp:163
void UpdateNode(const unsigned &rIndex, Node< SPACE_DIM > *pNode)
Definition: Element.cpp:89
void ResetIndex(unsigned index)
Definition: Element.cpp:104
c_vector< double, SPACE_DIM > CalculateXi(const ChastePoint< SPACE_DIM > &rTestPoint)
Definition: Element.cpp:290
void RegisterWithNodes()
Definition: Element.cpp:65
void MarkAsDeleted()
Definition: Element.cpp:74
void AddElement(unsigned index)
Definition: Node.cpp:269
Element(unsigned index, const std::vector< Node< SPACE_DIM > * > &rNodes, bool registerWithNodes=true)
Definition: Element.cpp:46
bool IncludesPoint(const ChastePoint< SPACE_DIM > &rTestPoint, bool strict=false)
Definition: Element.cpp:313
c_vector< double, SPACE_DIM+1 > CalculateInterpolationWeights(const ChastePoint< SPACE_DIM > &rTestPoint)
Definition: Element.cpp:228
c_vector< double, SPACE_DIM+1 > CalculateInterpolationWeightsWithProjection(const ChastePoint< SPACE_DIM > &rTestPoint)
Definition: Element.cpp:248
c_vector< double, 2 > CalculateMinMaxEdgeLengths()
Definition: Element.cpp:203