Real Coordinates of Gauss Points

[DanielB]

Dear all,

for my application I need to loop over all Gauss points of my surface model (build from 3D surface elements) and get their real coordinates. This I would need to do inside say a custom_utilities.h file, so outside a specific element or geometry description. Can somebody maybe suggest me a smart way to do that?

Already thanks in advance and best regards,
Daniel

[riccardo]

Hi Daniel

you should write a loop like the following

for(auto it = model_part.ElementsBegin(); it!=model_part.ElementsEnd(); it++)
{
   Geometry< Node<3> >& geom = it->GetGeometry();
   Matrix NContainer = geom.ShapeFunctionsValues(GeometryData::GI_GAUSS_2);

   
   for(unsigned int i=0; i<Ncontainer.size1(); i++)
   {
       array_1d<double,3> xgauss = ZeroVector(3);
       N = row(Ncontainer,i);

       for(unsigned int k=0; k<geom.size(); k++)
            xgauss += N*geom[k].Coordinates();

   }

}

[DanielB]

Hi Riccardo,

thanks a lot for this!!

In the meantime I figured out from a different place in the code, that I could also do it using preset functionalities of each geometry. So in particular I used:

for (ModelPart::ElementsContainerType::iterator elem_i = mr_model_part.ElementsBegin(); elem_i != mr_model_part.ElementsEnd(); ++elem_i)
        {
                const Element::GeometryType::IntegrationMethod integration_method = elem_i->GetIntegrationMethod();
           const Element::GeometryType::IntegrationPointsArrayType& integration_points = elem_i->GetGeometry().IntegrationPoints(integration_method);

           for ( unsigned int PointNumber = 0; PointNumber < integration_points.size(); PointNumber++ )
              NodeType::CoordinatesArrayType ip_coordinates = elem_i->GetGeometry().GlobalCoordinates(ip_coordinates, integration_points[PointNumber].Coordinates());
        }

Just let me know if you think this is a bad way of doing it.

Cheers,
Daniel

[DanielB]

A thing that popped up while trying both solutions: I have triangular elements / conditions where I need to loop over (each with a triangle_3d_3 geometry). When I call the function:

const Matrix& N_container = geom_i.ShapeFunctionsValues(integration_method);

I get an error telling me that the function "ShapeFunctionsValues" is not implemented for this particular triangle (so it calls the base class). For the quads however it works. So is there a particular reason for not including this function for triangular geometries? Or can you suggest maybe a way how I can do the above step such that it works for both triangles and quads?

In any case thanks a lot already and best wishes,
Daniel