65 assert(this->GetNumNodes() != 0);
67 if (this->GetNumNodes() <= 2)
84 for (
unsigned i=0; i<this->GetNumNodes(); i++)
86 mean_x += this->mNodes[i]->rGetLocation()[0];
87 mean_y += this->mNodes[i]->rGetLocation()[1];
89 mean_x /= this->GetNumNodes();
90 mean_y /= this->GetNumNodes();
92 double variance_x = 0;
93 double variance_y = 0;
94 double covariance_xy = 0;
96 for (
unsigned i=0; i<this->GetNumNodes(); i++)
98 variance_x += pow((this->mNodes[i]->rGetLocation()[0]-mean_x),2);
99 variance_y += pow((this->mNodes[i]->rGetLocation()[1]-mean_y),2);
100 covariance_xy += (this->mNodes[i]->rGetLocation()[0]-mean_x)*(this->mNodes[i]->rGetLocation()[1]-mean_y);
102 variance_x /= this->GetNumNodes();
103 variance_y /= this->GetNumNodes();
104 covariance_xy /= this->GetNumNodes();
107 double trace = variance_x+variance_y;
108 double det = variance_x*variance_y - covariance_xy*covariance_xy;
110 eig_max = 0.5*(trace+sqrt(trace*trace - 4*det));
111 eig_min = 0.5*(trace-sqrt(trace*trace - 4*det));
195 assert(eig_min >= 0);
196 assert(eig_max >= 0);
200 EXCEPTION(
"All nodes in an element lie in the same line/plane (2D/3D) so aspect ratio is infinite. This interferes with calculation of the Hamiltonian.");
203 return eig_max/eig_min;