Commit 16eadbc6 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

update coordinates 6.5x faster

parent 73d32f98
Pipeline #1293 passed with stage
in 26 minutes and 38 seconds
......@@ -34,7 +34,7 @@ static void _computeElementVector(const dgElementVector &ev, int integrationOrde
fullMatrix<double> coordEl;
for (int i = 0; i < ev.size(); ++i) {
ev.coordinates().getBlockProxy(i, coordEl);
elementInnerRadius(i) = fs.innerRadius(coordEl);
elementInnerRadius(i) = fs.innerRadius(coordEl, true);
fullMatrix<double> jac(3,3), ijac(3,3);
jac.setAll(0.);
fullMatrix<double> points(nbQP, 3);
......@@ -49,6 +49,7 @@ static void _computeElementVector(const dgElementVector &ev, int integrationOrde
st->applyJacobian(i, ev, xyz, uvw, jac, detJ(xi, i));
}
jac.invert(ijac);
for (int alpha = 0; alpha < fs.dimension; ++alpha) {
for (int a = 0; a < 3; a++) {
dXiDX (alpha, (i * nbQP + xi) * 3 + a) = ijac(a,alpha);
......@@ -104,21 +105,21 @@ static void _computeInterfaceVector(const dgInterfaceVector &iv, int integration
fullMatrix<double> jac(3,3), ijac(3,3);
int closureId = iv.closureId(i, iConn);
for (int xi = 0; xi < nbQP; xi++){
double coord[3] = {0, 0, 0};
double xyz[3] = {0, 0, 0};
for (int iPsi = 0; iPsi < nbPsi; iPsi++) {
for (int d = 0; d < fsElem.dimension; d++){
gReorder(fsElem.fullClosures[closureId][iPsi], d) = intMatrices.dPsi()(xi * fsElem.dimension + d, iPsi);
}
}
for (size_t iPsi = 0; iPsi < fsElem.closures[closureId].size(); iPsi++) {
for (int d = 0; d < fsElem.dimension; d++) {
coord[d] += intMatrices.psi()(xi, iPsi) * fsElem.points(fsElem.closures[closureId][iPsi], d);
xyz[d] += intMatrices.psi()(xi, iPsi) * coordEl(fsElem.closures[closureId][iPsi], d);
}
}
double j = fsElem.getJacobian(coordEl, gReorder, jac);
if (st) {
double coord[3] = {0, 0, 0};
double xyz[3] = {0, 0, 0};
for (size_t iPsi = 0; iPsi < fsElem.closures[closureId].size(); iPsi++) {
for (int d = 0; d < fsElem.dimension; d++) {
coord[d] += intMatrices.psi()(xi, iPsi) * fsElem.points(fsElem.closures[closureId][iPsi], d);
xyz[d] += intMatrices.psi()(xi, iPsi) * coordEl(fsElem.closures[closureId][iPsi], d);
}
}
st->applyJacobian(iv.elementId(i, iConn), iv.elementVector(iConn), xyz, coord, jac, j, closureId);
}
j = fabs(j);
......
......@@ -322,17 +322,38 @@ bool fullMatrix<double>::luSubstitute(const fullVector<double> &rhs, fullVector<
return false;
}
static double inv3x3(const fullMatrix<double> &m, fullMatrix<double> &i) {
i(0,0) = m(1,1)*m(2,2)-m(1,2)*m(2,1);
i(0,1) = m(2,1)*m(0,2)-m(2,2)*m(0,1);
i(0,2) = m(0,1)*m(1,2)-m(0,2)*m(1,1);
const double det = i(0,0)*m(0,0)+i(0,1)*m(1,0)+i(0,2)*m(2,0);
const double idet = 1./det;
i(0,0) *= idet;
i(0,1) *= idet;
i(0,2) *= idet;
i(1,0) = (m(2,0)*m(1,2)-m(1,0)*m(2,2))*idet;
i(1,1) = (m(0,0)*m(2,2)-m(2,0)*m(0,2))*idet;
i(1,2) = (m(1,0)*m(0,2)-m(0,0)*m(1,2))*idet;
i(2,0) = (m(1,0)*m(2,1)-m(1,1)*m(2,0))*idet;
i(2,1) = (m(2,0)*m(0,1)-m(2,1)*m(0,0))*idet;
i(2,2) = (m(0,0)*m(1,1)-m(0,1)*m(1,0))*idet;
return det;
}
template<>
bool fullMatrix<double>::invert(fullMatrix<double> &result) const
{
int M = size1(), N = size2(), lda = size1(), info;
int *ipiv = new int[std::min(M, N)];
if (result.size2() != M || result.size1() != N) {
if (result._own_data || !result._data)
result.resize(M,N,false);
else
Msg::Fatal("FullMatrix: Bad dimension, I cannot write in proxy");
}
if (M==3 && N ==3) {
return inv3x3(*this, result) != 0.;
}
int *ipiv = new int[std::min(M, N)];
result.setAll(*this);
F77NAME(dgetrf)(&M, &N, result._data, &lda, ipiv, &info);
if(info == 0){
......
......@@ -730,8 +730,19 @@ static void _planarQuad_xyz2xy(double *x, double *y, double *z, double *xn, doub
}
double nodalBasis::innerRadius(const fullMatrix<double> &x) const
double nodalBasis::innerRadius(const fullMatrix<double> &x, bool fastApproximate) const
{
if (fastApproximate && parentType != TYPE_LIN && parentType != TYPE_TRI && parentType != TYPE_TET) {
int nPrimaryNodes = primaryBasis().numNodes();
double inRadius = std::numeric_limits<double>::max();
for (int i = 0; i < nPrimaryNodes; ++i) {
for (int j = 0; j < i; ++j) {
const double vec[3] = {x(i,0)-x(j,0), x(i,1)-x(j,1),x(i,2)-x(j,2)};
inRadius = std::min(inRadius,vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
}
}
return sqrt(inRadius)/2;
}
switch(parentType) {
case TYPE_LIN :
{
......
......@@ -15,7 +15,7 @@ class nodalBasis {
int type, parentType, order, dimension, numFaces;
bool serendip;
fullMatrix<double> points;
double innerRadius(const fullMatrix<double> &nodes) const;
double innerRadius(const fullMatrix<double> &nodes, bool fastApproximate=false) const;
double volume(const fullMatrix<double> &nodes) const;
nodalBasis() {}
nodalBasis(int tag);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment