search

KratosMultiphysics / Kratos

Kratos Multiphysics (A.K.A Kratos) is a framework for building parallel multi-disciplinary simulation software. Modularity, extensibility and HPC are the main objectives. Kratos has BSD license and is written in C++ with extensive Python interface.
https://kratosmultiphysics.github.io/Kratos/
Other
1.03k stars 245 forks source link

Last calculated GP result is taken as result for all GPs in mesh #254

Closed peterjwilson closed 7 years ago

peterjwilson commented 7 years ago

Hi All,

Through benchmarking of my shell element stresses and forces (calculated on Gauss Points) I’ve found that certain meshes produce the problem that the last stress result calculated of the whole mesh is taken as the result for the whole mesh.

For example, the shell force XX result of a simply supported isotropic square plate with symmetry subject to a uniform transverse pressure.

Massimo’s existing thin triangle shell element: image

My new thick triangle shell element (with console output of the S_xx): image

The last calculated GP result in the console output is 0.267564, which is assigned to all elements in the mesh, despite the value varying as per the console output. Something strange is that my element doesn’t show this erroneous behaviour across other problems considered, for example a dome under self weight: image

Although it’s likely my element has some problems with it, it’s strange that the problem shows up in some meshes and not in others.

Has anyone experienced something like this, or has a hunch as to what could be going on here? Any help would be much appreciated.

jcotela commented 7 years ago

This is likely to be related to the definition of the default integration points (and/or quadrature rules) for your element. The first thing I would recommend would be to check the GetIntegrationMethod method for your element (I think this is the one used by the gid output) and check that it returns a quadrature with as many gauss points as you would like to print.

loumalouomega commented 7 years ago

Additionally I will recommend to get instead the Binary an Ascii file, this way you can see the defined coordinates of the GP and it values, some elements doesn't consider the standard quadrature and with the current API this can provide problems

peterjwilson commented 7 years ago

@jcotela I #define my integration method at the top of my cpp with the following line: #define OPT_INTEGRATION_METHOD Kratos::GeometryData::GI_GAUSS_1 Printing the GeometryType::IntegrationPointsArrayType from GetIntegrationMethod confirms I'm only using one gauss point. Is there a chance I'm using an outdated or incorrect method of setting my Gauss Point setup? I should note that I don't actually use any Kratos GP data while calculating quantities since it can all be analytically done at one point.

@loumalouomega Thanks for the tip, it provides some more clarity. The two relevant parts from the ASCII results are:

GaussPoints "tri1_element_gp" ElemType Triangle Number Of Gauss Points: 1 Natural Coordinates: Internal End gausspoints

and

Result "SHELL FORCE" "Kratos" 1.1 Matrix OnGaussPoints "tri1 element gp" ComponentNames "Sxx-SHELL FORCE", "Syy-SHELL FORCE", "Szz-SHELL FORCE", "Sxy-SHELL FORCE", "Syz-SHELL FORCE", "Sxz-SHELL FORCE" Values 1 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 2 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 3 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 4 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 5 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 6 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 7 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 8 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 9 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 10 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 11 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 12 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 13 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 14 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 15 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 16 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 17 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 18 0.267564416 0.206640214 0 0.126951098 424.213715 334.364685 End values

The first seems to look OK? But the second part confirms that the same result is written across all GPs.

peterjwilson commented 7 years ago

For those who might come across this problem in the future I solved it by adopting Massimo Petracca's approach in his ShellThinElement3D3N in the Structural Mechanics App. Namely setting the integration method to Kratos::GeometryData::GI_GAUSS_2 and then using the custom utilities Utilities::InterpToStandardGaussPoints(X) as defined in his element.