Hi,
whenever I try to use the Generate Quadrature Scheme Dictionary filter on my data with just 4 hexahedral finite elements having a different number of integration points each (8-node elements with full and reduced integration, 20-node elements with full and reduced integration), ParaView crashes. Here’s the vtk file:
elements.vtk (17.1 KB)
Is that because the Gauss points library doesn’t include integration point schemes of some of those elements ?
The error is:
Error: Cell type 25 found with no definition provided. Add a definition in /home/glow/dev/paraview/pv1/src/VTK/Filters/General/vtkQuadratureSchemeDictionaryGenerator.cxx. Aborting.
The segfault:
Thread 1 "paraview" received signal SIGSEGV, Segmentation fault.
vtkAOSDataArrayTemplate<long long>::GetTypedComponent (this=0x55556360e580, tupleIdx=1, comp=0) at /home/glow/dev/paraview/pv1/src/VTK/Common/Core/vtkAOSDataArrayTemplate.h:113
113 return this->Buffer->GetBuffer()[this->NumberOfComponents * tupleIdx + comp];
(gdb) bt
#0 vtkAOSDataArrayTemplate<long long>::GetTypedComponent (this=0x55556360e580, tupleIdx=1, comp=0) at /home/glow/dev/paraview/pv1/src/VTK/Common/Core/vtkAOSDataArrayTemplate.h:113
#1 0x00007fffedc8e33a in vtkGenericDataArray<vtkAOSDataArrayTemplate<long long>, long long>::SetTuple (this=0x555564b12d50, dstTupleIdx=0, srcTupleIdx=1, source=0x55556360e580) at /home/glow/dev/paraview/pv1/src/VTK/Common/Core/vtkGenericDataArray.txx:522
#2 0x00007fffedc8b43c in vtkAOSDataArrayTemplate<long long>::SetTuple (this=0x555564b12d50, dstTupleIdx=0, srcTupleIdx=1, source=0x55556360e580) at /home/glow/dev/paraview/pv1/src/VTK/Common/Core/vtkAOSDataArrayTemplate.h:201
#3 0x00007fffedc8e43d in vtkGenericDataArray<vtkAOSDataArrayTemplate<long long>, long long>::InsertTuple (this=0x555564b12d50, i=0, j=1, source=0x55556360e580) at /home/glow/dev/paraview/pv1/src/VTK/Common/Core/vtkGenericDataArray.txx:679
#4 0x00007fffedc8b698 in vtkAOSDataArrayTemplate<long long>::InsertTuple (this=0x555564b12d50, dstTupleIdx=0, srcTupleIdx=1, source=0x55556360e580) at /home/glow/dev/paraview/pv1/src/VTK/Common/Core/vtkAOSDataArrayTemplate.h:209
#5 0x00007fffef4afa58 in vtkDataSetAttributes::CopyTuple (this=0x555564b093f0, fromData=0x55556360e580, toData=0x555564b12d50, fromId=1, toId=0) at /home/glow/dev/paraview/pv1/src/VTK/Common/DataModel/vtkDataSetAttributes.cxx:1172
#6 0x00007fffef4ae7c8 in vtkDataSetAttributes::CopyData (this=0x555564b093f0, fromPd=0x555564b0c1b0, fromId=1, toId=0) at /home/glow/dev/paraview/pv1/src/VTK/Common/DataModel/vtkDataSetAttributes.cxx:924
#7 0x00007fffdf0aafd4 in vtkUnstructuredGridGeometryFilter::RequestData (this=0x555563c10470, inputVector=0x555564b206f0, outputVector=0x555564b20610) at /home/glow/dev/paraview/pv1/src/VTK/Filters/Geometry/vtkUnstructuredGridGeometryFilter.cxx:1270
#8 0x00007ffff040876e in vtkUnstructuredGridBaseAlgorithm::ProcessRequest (this=0x555563c10470, request=0x555564b1b000, inputVector=0x555564b206f0, outputVector=0x555564b20610) at /home/glow/dev/paraview/pv1/src/VTK/Common/ExecutionModel/vtkUnstructuredGridBaseAlgorithm.cxx:69
#9 0x00007ffff028be4b in vtkExecutive::CallAlgorithm (this=0x555564a8c300, request=0x555564b1b000, direction=1, inInfo=0x555564b206f0, outInfo=0x555564b20610) at /home/glow/dev/paraview/pv1/src/VTK/Common/ExecutionModel/vtkExecutive.cxx:735
#10 0x00007ffff02821d8 in vtkDemandDrivenPipeline::ExecuteData (this=0x555564a8c300, request=0x555564b1b000, inInfo=0x555564b206f0, outInfo=0x555564b20610) at /home/glow/dev/paraview/pv1/src/VTK/Common/ExecutionModel/vtkDemandDrivenPipeline.cxx:462
#11 0x00007ffff0275411 in vtkCompositeDataPipeline::ExecuteData (this=0x555564a8c300, request=0x555564b1b000, inInfoVec=0x555564b206f0, outInfoVec=0x555564b20610) at /home/glow/dev/paraview/pv1/src/VTK/Common/ExecutionModel/vtkCompositeDataPipeline.cxx:163
#12 0x00007ffff0281785 in vtkDemandDrivenPipeline::ProcessRequest (this=0x555564a8c300, request=0x555564b1b000, inInfoVec=0x555564b206f0, outInfoVec=0x555564b20610) at /home/glow/dev/paraview/pv1/src/VTK/Common/ExecutionModel/vtkDemandDrivenPipeline.cxx:261
#13 0x00007ffff03d8032 in vtkStreamingDemandDrivenPipeline::ProcessRequest (this=0x555564a8c300, request=0x555564b1b000, inInfoVec=0x555564b206f0, outInfoVec=0x555564b20610) at /home/glow/dev/paraview/pv1/src/VTK/Common/ExecutionModel/vtkStreamingDemandDrivenPipeline.cxx:344
#14 0x00007ffff0281fbd in vtkDemandDrivenPipeline::UpdateData (this=0x555564a8c300, outputPort=0) at /home/glow/dev/paraview/pv1/src/VTK/Common/ExecutionModel/vtkDemandDrivenPipeline.cxx:419
This filter does not support Quadratic Hexahedron, but it shouldn’t segfault either. Please open an issue: https://gitlab.kitware.com/paraview/paraview/-/issues
A possible workaround is to Tetrahadrlize your dataset, but I’m not sure that is what you want.
FYI @jfausty
Sure, I’ll do it.
Would it be difficult to add it ? Maybe I could even help. What are the currently included element types ? I wanted to test linear hexahedron with full and reduced integration, quadratic hex with full and reduced integration and then also linear and quadratic tetrahedron. 2D elements are excluded from those tests for now because CalculiX expands them to solids and uses some non-standard integration point scheme. Also, how are their schemes defined in that dictionary (I guess that this file is the one to expand: VTK/Filters/General/vtkQuadratureSchemeDictionaryGenerator.cxx at master · Kitware/VTK · GitHub) ?
It’s right here: https://gitlab.kitware.com/vtk/vtk/-/blob/master/Filters/General/vtkQuadratureSchemeDictionaryGenerator.cxx#L257
Please take a look here to learn how to contribute to VTK: https://gitlab.kitware.com/vtk/vtk/-/blob/master/Documentation/dev/git/develop.md#quick-start-guide
Let me know if you need any help.
FYI @paraview
I see, so only linear and quadratic elements of the following types are supported: triangle, quadrilateral, tetrahedron and thus hexahedrons are not supported at all. Would it be possible to also add support for reduced integration versions of some of those elements - is it possible to determine whether full or reduced integration scheme is used ? Or maybe a user-activated option for that would suffice.
Thanks, I’ll take a look at this. Do you know how those definitions for each element type are made ? For example, for linear quad it looks like this:
double W_Q_42_A[] = { 6.22008467928145e-01, 1.66666666666667e-01, 4.46581987385206e-02,
1.66666666666667e-01, 1.66666666666667e-01, 4.46581987385206e-02, 1.66666666666667e-01,
6.22008467928145e-01, 1.66666666666667e-01, 6.22008467928145e-01, 1.66666666666667e-01,
4.46581987385206e-02, 4.46581987385206e-02, 1.66666666666667e-01, 6.22008467928145e-01,
1.66666666666667e-01 };
Those numbers are described as “shape functions weights” but I wonder if they can be found in this form in the literature. If yes then it shouldn’t be difficult to add support for other element types.
Hi @FEngineer,
Unfortunately, I can’t seem to trace back the original definition of these arrays. I imagine they represent the matrices used to interpolate values from Lagrangian nodes to quadrature points (which do not describe the rules completely since we are still missing the weights).
The W_Q_42_A one seems to have 4x4 entries which would mean it has 4 quadrature points and, assuming it is a tensor product quadrature, should be exact up to 3rd order polynomials. I don’t know whether this corresponds to your definition of “full” vs. “reduced” integration. For general orhtotope (line, quad, hexahedra) finite elements, you usually have to use quadrature rules with orders of 4p to be exact (where p is the polynomial order of your element).
Additionally, reading the documentation from https://gitlab.kitware.com/vtk/vtk/-/blob/master/Filters/General/vtkQuadratureSchemeDictionaryGenerator.h, it would seem that this class was destined to be used as a test utility and should not be used for application development.
May I ask in what context are these quadrature rules of use to you in ParaView?
Hope that helps!
Best regards,
Julien
Thank you for the reply.
That would make sense but I’m afraid that it will be hard to find such matrices in the literature. For some reason, they use very specific numbers like 6.22008467928145e-01.
Right, but I wonder why e.g. quadratic quad (W_QQ_93_A) is described as 8x9. It has 9 integration points.
Reduced integration:
That’s strange, this functionality seems to be working well and is very useful. Maybe that note is outdated.
They are very useful for the postprocessing of FEA results in which stresses are given as interpolated and averaged at nodes (which is often the case). For the most accurate results, we check the values at the integration points directly.
Hi @FEngineer,
Glad I could help.
Right, but I wonder why e.g. quadratic quad (
W_QQ_93_A) is described as 8x9. It has 9 integration points.
That is correct, the quadratic quad has 8 nodes (4 on vertices and 4 on the midpoints of its edges) and the 9 integration points correspond to a quadrature rule with a maximal polynomial order of 5.
Reduced integration:
linear quad - 1 int point (instead of 4)
quadratic quad - 4 int points (instead of 9)
If you were to improve this filter by contributing, the reduced quadrature rules would need to be controlled by a a new filter property in order to resolve the ambiguity (for example, if you have a quadratic quad in your data set it would use the 9 point quadrature rule by default now).
They are very useful for the postprocessing of FEA results in which stresses are given as interpolated and averaged at nodes (which is often the case). For the most accurate results, we check the values at the integration points directly.
I see, in that case the way forward seems to be indeed adding a contribution to VTK (following https://gitlab.kitware.com/vtk/vtk/-/blob/master/Documentation/dev/git/develop.md#quick-start-guide) and then update the VTK changes in ParaView master. If you are looking to generate your meshes with quadrature points included, for best results, I would advise you to use the same quadrature rules as defined in your FE software. If that is not possible, most of the time, quadrature rules for quadrilaterals and hexahedra are tensor products of the classic 1D Gaussian quadrature rules, which can be generated in a straightforward manner.
Best regards,
Julien
Now it makes sense, I forgot that the first number corresponds to the number of nodes.
Yes, I was thinking about adding a switch that would allow a user to change from the default full integration to reduced integration. But support for hex elements should be added first.
I will try to do that. The biggest problem now is that I don’t know where to find matrices defining integration schemes in the same way as it’s done in the vtkQuadratureSchemeDictionaryGenerator.cxx file.
I will try to do that. The biggest problem now is that I don’t know where to find matrices defining integration schemes in the same way as it’s done in the vtkQuadratureSchemeDictionaryGenerator.cxx file.
If there are no constraints on your side, the simplest way I can see is to take the set of Gaussian quadrature points on the [-1, 1] interval corresponding to your polynomial order, let us denote them
\mathcal{G}^{1}_{N} = \{p_{i}, \quad i = 1,\ldots,N\}
and run the tensor product to generate higher dimensional points
\mathcal{G}^{d}_{N} = \{(p_{i}, p_{j}, \ldots, p_{k}) \in (\mathcal{G}^{1}_{N})^d, \quad (i, j, \ldots, k) = \{1, \ldots, N\}^d \}
Once you have the points, the interpolation matrix is computed by interpolating each one of your basis functions on each one of the points to get a (N^d, M) matrix where M is the number of basis functions in your element (usually the number of nodes).
@jfausty I’ll try to figure it out. But there’s one more thing bothering me - why is a linear triangular element defined as having 3 integration points (3x3) ? Shouldn’t it have 1 such point or is it just a matter of different formulations ?
It is a matter of formulations. A single integration point would approximate values inside the cell as a constant. As the number of integration points increases, the order of the polynomial used to approximate the finite element solution inside the cell increases. Typically, higher order polynomials limit the error bounds on the integral. Common integration point placements are Gauss or Gauss-Lobatto points that minimize the error under different conditions.
So you can integrate linear functions using a single integration point and be exact. The reason FE software uses higher than the cell order quadrature schemes is that in Galerkin weak formulations of problems (and their subsequent FE discretizations) you end up with terms that look like:
\int_{\Omega}\phi_{i}\phi_{j}d\Omega
where \{\phi_{i}\} are your basis functions. As such, what needs to integrated is no longer of order p, the order of your basis functions, but of order 2p necessitating a quadrature rule that can integrate at least 2p polynomials in your cell. For quads and hexahedra, the actual order of the basis functions of a cell of order p can rise up to 2p because of the cross coordinate terms and thus necessitate quadrature rules that can account for at least 4p.
@jfausty I’ll try to figure it out. But there’s one more thing bothering me - why is a linear triangular element defined as having 3 integration points (3x3) ? Shouldn’t it have 1 such point or is it just a matter of different formulations ?
I believe this choice was made to be consistent with FE software which actually use this order of quadrature rule on linear cells for the above mentioned reasons. The three point quadrature rule in the triangle can integrate quadratic polynomials.
Thank you for the replies. In such a case, an additional option to switch from a linear triangle with 3 int points to its version with 1 int point (used by Abaqus and CalculiX) could be useful. Same for linear tetrahedrons since ParaView uses 4 integration points for those elements while Abaqus and CalculiX use only 1.
I think that I am close to figuring out how those matrices should look for hex elements but I would appreciate your additional help. I’ve found similar file in CalculiX source code: CalculiX/gauss.f at master · Dhondtguido/CalculiX · GitHub
There you can find two matrices for each element type, they are named gauss (integration point locations) and weight (weight functions). The problem is that I don’t know how those two matrices can be used together to obtain a matrix like the one in ParaView. It can’t be just a simple product of those two matrices.
For example, let’s focus on quadratic tetrahedron since it’s the same in CalculiX and ParaView. Matrix for this element in ParaView is:
{ 1.56250000000000e-01, -9.37500000000000e-02, -9.37500000000000e-02,
-9.37500000000000e-02, 3.12500000000000e-01, 6.25000000000000e-02, 3.12500000000000e-01,
3.12500000000000e-01, 6.25000000000000e-02, 6.25000000000000e-02, -9.37500000000000e-02,
7.03125000000000e-02, -1.17187500000000e-01, -9.37500000000000e-02, 2.81250000000000e-01,
4.21875000000000e-01, 9.37500000000000e-02, 6.25000000000000e-02, 2.81250000000000e-01,
9.37500000000000e-02, -9.37500000000000e-02, -1.17187500000000e-01, 7.03125000000000e-02,
-9.37500000000000e-02, 9.37500000000000e-02, 4.21875000000000e-01, 2.81250000000000e-01,
6.25000000000000e-02, 9.37500000000000e-02, 2.81250000000000e-01, -9.37500000000000e-02,
-5.46875000000000e-02, -5.46875000000000e-02, 3.75000000000000e-01, 3.12500000000000e-02,
1.56250000000000e-02, 3.12500000000000e-02, 3.75000000000000e-01, 1.87500000000000e-01,
1.87500000000000e-01 }
while matrices in CalculiX (using Fortran) are:
reshape((/
& 0.138196601125011d0,0.138196601125011d0,0.138196601125011d0,
& 0.585410196624968d0,0.138196601125011d0,0.138196601125011d0,
& 0.138196601125011d0,0.585410196624968d0,0.138196601125011d0,
& 0.138196601125011d0,0.138196601125011d0,0.585410196624968d0/),
& (/3,4/))
(/
& 0.041666666666667d0,0.041666666666667d0,0.041666666666667d0,
& 0.041666666666667d0/)
Also, here’s part of the book used as a theory manual for CalculiX:
Since e.g. (5-sqrt(5))/20=0.1381966 and (5+(3*sqrt(5)))/20=0.5854102 and 1/24=0.0416667, the numbers in both matrices from CalculiX source code make sense. But I still don’t know how to get the matrix like the one in ParaView from those two matrices from CalculiX (especially since they don’t account for the number of nodes). Any ideas ? Knowing how it works for tetrahedral elements, I could do the same for hex elements.
Hi @FEngineer,
Great!
Indeed, in CalculiX it seems that the quadrature rules are in a more conventional format.
From what I understood from the quadrature definitions in VTK, the matrices inputted into the vtkQuadratureSchemeDefinition object are the values of the shape functions at the intergration points directly. Let’s note this matrix for a cell M_{interp} with size (N_{integrationpoints}, N_{numbernodes}).
In the VTK context, the coordinates of the integration points can be calculated from these matrices by matrix product with the coordinates of the nodes in matrix form C_{nodes} of shape (N_{numbernodes}, d) where d is the topological dimension of the cell. As such:
C_{integrationpoints} = M_{interp}C_{nodes}
similarly to any other field defined in the cell.
As such, to retrieve this M_{interp} matrix from your CalculiX cell, you have but to evaluate you shape functions at the integration points as so:
M_{interp} = \begin{pmatrix} \phi^{shape}_{0}(X_{integration}^{0}) & \ldots & \phi^{shape}_{N_{numbernodes}}(X^{0}_{integration}) \\ & \ldots & \\ \phi^{shape}_{0}(X_{integration}^{N_{integrationpoints}}) & \ldots & \phi^{shape}_{N_{numbernodes}}(X^{N_{integrationpoints}}_{integration}) \end{pmatrix}
The weight information is not included here as the weight information at each of the points is not included in these matrices. The quadrature definition is thus incomplete but sufficient to transport field data from the nodes of your element to its integration points. If you are looking to also integrate those fields, you will have to include the weights as well which the vtkQuadratureSchemeDefinition supports.
Thank you.
As I understand, X_{integration} means locations of the integration points so basically gauss3d5 matrix from CalculiX, right ? And it should be multiplied by shape functions \phi^{shape} (unfortunately, not given in that CalculiX file) so essentially functions of local coordinates of integration points in this case if I understand it correctly. I just don’t get the matrix sizes. According to the book, there are 4 coordinates for each of the 4 integration points but the gauss3d5 matrix in CalculiX is of 3x4 size for some reason. And there are 10 shape functions (since quadratic tetrahedron has 10 nodes). The resulting matrix needs to be of (4,10) size like W_QE_42_A in ParaView. So how would the shape function matrix look like - should it be 4x4 somehow ?
So the W_QE_42_A matrix in vtkQuadratureSchemeDictionaryGenerator.cxx is M_{interp} but apart from that weights need to be included in an additional file named vtkQuadratureSchemeDefinition.cxx ?
As I understand, X_{integration} means locations of the integration points so basically gauss3d5 matrix from CalculiX, right ?
That’s right
And it should be multiplied by shape functions \phi^{shape} (unfortunately, not given in that CalculiX file) so essentially functions of local coordinates of integration points in this case if I understand it correctly.
That’s where I should have maybe explained it better. The M_{interp} matrix components are the values of the shape functions at the integration points (this is not a matrix multiplication operation). This is a commonly precomputed matrix in FE calculations that is used in the assembly of the local system matrices. A CalculiX cell should have a method or a member which depends on the quadrature rule and computes this matrix.
This subject is more related to the internals of CalculiX and not ParaView though so I am afraid I won’t have much more information for you on that. Do you have any way to contact their support or maintainers?
I just don’t get the matrix sizes. According to the book, there are 4 coordinates for each of the 4 integration points but the gauss3d5 matrix in CalculiX is of 3x4 size for some reason.
I believe the book is using barycentric coordinates to describe the location while the text file has Cartesian coordinates. I could be wrong though. Another question for the maintainers of the CalculiX code.
So the W_QE_42_A matrix in vtkQuadratureSchemeDictionaryGenerator.cxx is Minterp but apart from that weights need to be included in an additional file named vtkQuadratureSchemeDefinition.cxx ?
The vtkQuadratureSchemeDictionaryGenerator generates a dictionary of objects of type vtkQuadratureSchemeDefinition. The vtkQuadratureSchemeDefinition object has methods for setting the quadrature weights in addition to the shape function matrix (https://gitlab.kitware.com/vtk/vtk/-/blob/master/Common/DataModel/vtkQuadratureSchemeDefinition.h). However, the only reason you would need to set the weights was if you were performing integrals from within ParaView. This does not seem part of your goals. As such, I do not believe you need to include the weight information in your development.
As this discussion is veering away from the strictly VTK/ParaView side of things I would advise you to try to contact the maintainers of the CalculiX code for more detailed explanations of their cell structures.
Thanks again. I will take a closer look at CalculiX source code and maybe ask someone where to find the necessary objects.
I should have guessed by those parentheses.
To sum up:
I need a 4x10 matrix containing shape function values at the integration points for the quadratic tetrahedron element. This matrix should be equivalent to the W_QE_42_A one in the vtkQuadratureSchemeDictionaryGenerator.cxx file and I don’t have to worry about weights.
