Commit fe3ff3b2 authored by Berk Geveci's avatar Berk Geveci

Fixed quadratic quad and biquadratic quad derivatives.

Fixed bug #13455:
vtkCellDerivatives computes wrong results on quadratic quadrilateral
cells.

Adapted fix from Simon Triebenbacher. Thanks Simon.

Change-Id: I49d0064adcdadd5780e24a6a10e039159a3dc850
parent 562e93f9
......@@ -13,5 +13,6 @@ vtk_add_test_python(
TestConvexPointSet.py
otherDataSetAttributes.py,NO_DATA,NO_VALID,NO_RT
quadricCut.py,NO_RT
quadraticQuadDeriv.py,NO_VALID,NO_RT
TestICPTransform.py,NO_RT
)
#!/usr/bin/env python
import vtk
from vtk.test import Testing
from vtk.util.misc import vtkGetDataRoot
VTK_DATA_ROOT = vtkGetDataRoot()
from vtk.test import Testing
class TestCommand(Testing.vtkTest):
def _test(self, fname):
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(VTK_DATA_ROOT + fname)
elev = vtk.vtkElevationFilter()
elev.SetInputConnection(reader.GetOutputPort())
elev.SetLowPoint(-0.05, 0.05, 0)
elev.SetHighPoint(0.05, 0.05, 0)
grad = vtk.vtkGradientFilter()
grad.SetInputConnection(elev.GetOutputPort())
grad.SetInputArrayToProcess(0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "Elevation")
grad.Update()
vals = (10, 0, 0)
for i in range(3):
r = grad.GetOutput().GetPointData().GetArray("Gradients").GetRange(i)
self.assertTrue(abs(r[0] - vals[i]) < 1E-4)
self.assertTrue(abs(r[1] - vals[i]) < 1E-4)
elev.SetLowPoint(0.05, -0.05, 0)
elev.SetHighPoint(0.05, 0.05, 0)
grad.Update()
vals = (0, 10, 0)
for i in range(3):
r = grad.GetOutput().GetPointData().GetArray("Gradients").GetRange(i)
self.assertTrue(abs(r[0] - vals[i]) < 1E-4)
self.assertTrue(abs(r[1] - vals[i]) < 1E-4)
def testQuadraticQuad(self):
self._test("/Data/Disc_QuadraticQuads_0_0.vtu")
def testBiQuadraticQuad(self):
self._test("/Data/Disc_BiQuadraticQuads_0_0.vtu")
if __name__ == "__main__":
Testing.main([(TestCommand, 'test')])
......@@ -335,34 +335,73 @@ vtkBiQuadraticQuad::Derivatives (int vtkNotUsed (subId),
double pcoords[3], double *values,
int dim, double *derivs)
{
double x0[3], x1[3], x2[3], deltaX[3], weights[9];
int i, j;
double sum[3], weights[9];
double functionDerivs[18];
double elemNodes[9][3];
double *J[3], J0[3], J1[3], J2[3];
double *JI[3], JI0[3], JI1[3], JI2[3];
this->Points->GetPoint (0, x0);
this->Points->GetPoint (1, x1);
this->Points->GetPoint (2, x2);
for(int i = 0; i<9; i++)
{
this->Points->GetPoint(i, elemNodes[i]);
}
this->InterpolationFunctions (pcoords, weights);
this->InterpolationDerivs (pcoords, functionDerivs);
this->InterpolationFunctions(pcoords,weights);
this->InterpolationDerivs(pcoords,functionDerivs);
for (i = 0; i < 3; i++)
// Compute transposed Jacobian and inverse Jacobian
J[0] = J0; J[1] = J1; J[2] = J2;
JI[0] = JI0; JI[1] = JI1; JI[2] = JI2;
for(int k = 0; k<3; k++)
{
deltaX[i] = x1[i] - x0[i] - x2[i];
J0[k] = J1[k] = 0.0;
}
for (i = 0; i < dim; i++)
for(int i = 0; i<9; i++)
{
for (j = 0; j < 3; j++)
for(int j = 0; j<2; j++)
{
if (deltaX[j] != 0)
for(int k = 0; k<3; k++)
{
derivs[3 * i + j] = (values[2 * i + 1] - values[2 * i]) / deltaX[j];
J[j][k] += elemNodes[i][k] * functionDerivs[j*9+i];
}
else
}
}
// Compute third row vector in transposed Jacobian and normalize it, so that Jacobian determinant stays the same.
vtkMath::Cross(J0,J1,J2);
if ( vtkMath::Normalize(J2) == 0.0 || !vtkMath::InvertMatrix(J,JI,3)) //degenerate
{
for (int j=0; j < dim; j++ )
{
for (int i=0; i < 3; i++ )
{
derivs[3 * i + j] = 0;
derivs[j*dim + i] = 0.0;
}
}
return;
}
// Loop over "dim" derivative values. For each set of values,
// compute derivatives
// in local system and then transform into modelling system.
// First compute derivatives in local x'-y' coordinate system
for (int j=0; j < dim; j++ )
{
sum[0] = sum[1] = sum[2] = 0.0;
for (int i=0; i < 9; i++) //loop over interp. function derivatives
{
sum[0] += functionDerivs[i] * values[dim*i + j];
sum[1] += functionDerivs[9 + i] * values[dim*i + j];
}
// dBydx = sum[0]*JI[0][0] + sum[1]*JI[0][1];
// dBydy = sum[0]*JI[1][0] + sum[1]*JI[1][1];
// Transform into global system (dot product with global axes)
derivs[3*j] = sum[0]*JI[0][0] + sum[1]*JI[0][1];
derivs[3*j + 1] = sum[0]*JI[1][0] + sum[1]*JI[1][1];
derivs[3*j + 2] = sum[0]*JI[2][0] + sum[1]*JI[2][1];
}
}
......
......@@ -452,34 +452,73 @@ void vtkQuadraticQuad::Derivatives(int vtkNotUsed(subId),
double pcoords[3], double *values,
int dim, double *derivs)
{
double x0[3], x1[3], x2[3], deltaX[3], weights[8];
int i, j;
double sum[3], weights[8];
double functionDerivs[16];
double elemNodes[8][3];
double *J[3], J0[3], J1[3], J2[3];
double *JI[3], JI0[3], JI1[3], JI2[3];
this->Points->GetPoint(0, x0);
this->Points->GetPoint(1, x1);
this->Points->GetPoint(2, x2);
for(int i = 0; i<8; i++)
{
this->Points->GetPoint(i, elemNodes[i]);
}
this->InterpolationFunctions(pcoords,weights);
this->InterpolationDerivs(pcoords,functionDerivs);
for (i=0; i<3; i++)
// Compute transposed Jacobian and inverse Jacobian
J[0] = J0; J[1] = J1; J[2] = J2;
JI[0] = JI0; JI[1] = JI1; JI[2] = JI2;
for(int k = 0; k<3; k++)
{
deltaX[i] = x1[i] - x0[i] - x2[i];
J0[k] = J1[k] = 0.0;
}
for (i=0; i<dim; i++)
for(int i = 0; i<8; i++)
{
for (j=0; j<3; j++)
for(int j = 0; j<2; j++)
{
if ( deltaX[j] != 0 )
for(int k = 0; k<3; k++)
{
derivs[3*i+j] = (values[2*i+1] - values[2*i]) / deltaX[j];
J[j][k] += elemNodes[i][k] * functionDerivs[j*8+i];
}
else
}
}
// Compute third row vector in transposed Jacobian and normalize it, so that Jacobian determinant stays the same.
vtkMath::Cross(J0,J1,J2);
if ( vtkMath::Normalize(J2) == 0.0 || !vtkMath::InvertMatrix(J,JI,3)) //degenerate
{
for (int j=0; j < dim; j++ )
{
for (int i=0; i < 3; i++ )
{
derivs[3*i+j] = 0;
derivs[j*dim + i] = 0.0;
}
}
return;
}
// Loop over "dim" derivative values. For each set of values,
// compute derivatives
// in local system and then transform into modelling system.
// First compute derivatives in local x'-y' coordinate system
for (int j=0; j < dim; j++ )
{
sum[0] = sum[1] = sum[2] = 0.0;
for (int i=0; i < 8; i++) //loop over interp. function derivatives
{
sum[0] += functionDerivs[i] * values[dim*i + j];
sum[1] += functionDerivs[8 + i] * values[dim*i + j];
}
// dBydx = sum[0]*JI[0][0] + sum[1]*JI[0][1];
// dBydy = sum[0]*JI[1][0] + sum[1]*JI[1][1];
// Transform into global system (dot product with global axes)
derivs[3*j] = sum[0]*JI[0][0] + sum[1]*JI[0][1];
derivs[3*j + 1] = sum[0]*JI[1][0] + sum[1]*JI[1][1];
derivs[3*j + 2] = sum[0]*JI[2][0] + sum[1]*JI[2][1];
}
}
......
5fb60892fb8d76bd5de9bd1b9f4003e3
e9e7cd584e083a5b7c807eb0accbfad6
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