[Core] Line/triangle IsInside removing projection

[philbucher]

Description
This PR removes the projection that is happening in the IsInside function of Line2D2 and Triangle3D3
Reasons:

• Inconsistent with other geometries, there only the basic check is done
• Projection is sort of an informative/safety feature, but it punishes users that use the function correctly as then projection is unnecessary!
• Affects the shared memory parallelism very much. Even for my simple example as shown the performance takes a serious hit!!! I don't even want to think what can happen in bigger examples. I think this is manly related to the Logger macros that have a shared counter (IIRC @roigcarlo)

results for 1e7 executions with the projection

	Line (time [s])	Triangle (time [s])
serial	0.251562	0.709768
parallel (2 threads)	0.432373	0.578865
parallel (3 threads)	0.333716	0.453672
parallel (4 threads)	0.323167	0.39459

results for 1e7 executions without the projection (here it pretty much scales linear, as expected)

	Line (time [s])	Triangle (time [s])
serial	0.134439	0.424116
parallel (2 threads)	0.0565617	0.212631
parallel (3 threads)	0.0381214	0.144309
parallel (4 threads)	0.0295059	0.110287

[philbucher]

for reference the code I used to profile:

#include "testing/testing.h"
#include "utilities/builtin_timer.h"
#include "utilities/openmp_utils.h"
#include "geometries/line_2d_2.h"
#include "geometries/triangle_3d_3.h"

namespace Kratos {
namespace Testing {
using GeometryType = Geometry<Point>;
using GeometryPointerType = Kratos::shared_ptr<GeometryType>;
using CoordinatesArrayType = Geometry<Point>::CoordinatesArrayType;
static constexpr std::size_t num_evals = 1e7;

void RunPerformanceTest(
    const std::string& rTestName,
    GeometryPointerType pGeom,
    const CoordinatesArrayType& rPointInside,
    const CoordinatesArrayType& rPointOutside)
{
    CoordinatesArrayType result;

    BuiltinTimer timer_i;
    for (std::size_t i=0; i<num_evals; ++i) {
        pGeom->IsInside(rPointInside, result);
    }
    const double time_i = timer_i.ElapsedSeconds();

    BuiltinTimer timer_i_p;
    #pragma omp parallel for firstprivate(rPointInside), private(result)
    for (std::size_t i=0; i<num_evals; ++i) {
        pGeom->IsInside(rPointInside, result);
    }
    const double time_i_p = timer_i_p.ElapsedSeconds();

    BuiltinTimer timer_o;
    for (std::size_t i=0; i<num_evals; ++i) {
        pGeom->IsInside(rPointOutside, result);
    }
    const double time_o = timer_o.ElapsedSeconds();

    BuiltinTimer timer_o_p;
    #pragma omp parallel for firstprivate(rPointOutside), private(result)
    for (std::size_t i=0; i<num_evals; ++i) {
        pGeom->IsInside(rPointOutside, result);
    }
    const double time_o_p = timer_o_p.ElapsedSeconds();

    std::cout << "\n\nRunning Test for: " << rTestName << " with " << OpenMPUtils::GetNumThreads() << " threads\n"
              << "With " << num_evals << " number of evaluations the results are:\n"
              << "Inside, serial    : " << time_i << " [s]\n"
              << "Inside, parallel  : " << time_i_p << " [s]\n"
              << "Outside, serial   : " << time_o << " [s]\n"
              << "Outside, parallel : " << time_o_p << " [s]\n\n" << std::endl;

}

KRATOS_TEST_CASE_IN_SUITE(Line2D2, KratosCoreFastSuite)
{
    GeometryPointerType geom = Kratos::make_shared<Line2D2<Point>>(
        Kratos::make_shared<Point>(0.0, 0.0, 0.0),
        Kratos::make_shared<Point>(1.0, 0.0, 0.0)
        );

    CoordinatesArrayType point_inside;
    point_inside[0] = 0.5;
    point_inside[1] = 0.0;
    point_inside[2] = 0.0;

    CoordinatesArrayType point_outside;
    point_outside[0] = 0.5;
    point_outside[1] = 0.1;
    point_outside[2] = 0.0;
    CoordinatesArrayType result;

    RunPerformanceTest("Line2D2", geom, point_inside, point_outside);
}

KRATOS_TEST_CASE_IN_SUITE(Triangle3D3, KratosCoreFastSuite)
{
    GeometryPointerType geom = Kratos::make_shared<Triangle3D3<Point>>(
        Kratos::make_shared<Point>(0.0, 0.0, 0.0),
        Kratos::make_shared<Point>(1.0, 0.0, 0.0),
        Kratos::make_shared<Point>(1.0, 1.0, 0.0)
        );

    CoordinatesArrayType point_inside;
    point_inside[0] = 0.1;
    point_inside[1] = 0.1;
    point_inside[2] = 0.0;

    CoordinatesArrayType point_outside;
    point_outside[0] = 0.1;
    point_outside[1] = 0.1;
    point_outside[2] = 0.1;
    CoordinatesArrayType result;

    RunPerformanceTest("Triangle3D3", geom, point_inside, point_outside);
}

} // namespace Testing
}  // namespace Kratos.

[loumalouomega]

Blocking as I need to check if it breaks Contact (I modified the code after we included the projection)

[roigcarlo]

@philbucher Indeed that macro causes an important performance hit if combined with openmp. I will prepare a PR changing the rest of the uses in Kratos, as they may have similar problems.

Also probably is good idea to only enable macros with counters while in fulldebug

[pooyan-dadvand]

@KratosMultiphysics/technical-committee agrees with this change due to the following reasons:

• It is more coherent with the other geometries
• If the point is outside and one wanted to check the projected (using ProjectionPoint method), it can be done by first projecting the point and then check the IsInside
• In the current implementation one cannot check if point is simply outside because it is out of plane.

@loumalouomega we would wait until you verify and confirm that it works for you.

[loumalouomega]

I already checked the contact, but it breaks the Meshing, see tests

[philbucher]

I already checked the contact, but it breaks the Meshing, see tests

yes we saw that but can you check if you can refactor?
I.e. doing what was done in the IsInside function in the code before calling IsInside

[loumalouomega]

I already checked the contact, but it breaks the Meshing, see tests

yes we saw that but can you check if you can refactor?
I.e. doing what was done in the IsInside function in the code before calling IsInside

I will

[loumalouomega]

There are fails, I need to do some changes, sorry

[loumalouomega]

I have not finished with the changes

[tteschemacher]

In a recent technical committee meeting we agreed to have the Is inside as a function including the projection. If the projection is done beforehand it was suggested to use IsInsideLocalSpace wouldn't it be the most consistent moving towards this interface.

[philbucher]

In a recent technical committee meeting we agreed to have the Is inside as a function including the projection. If the projection is done beforehand it was suggested to use IsInsideLocalSpace wouldn't it be the most consistent moving towards this interface.

Hm this sounds like a lot of work and is not applicable for e.g. solids
Now thinking about it the other approach seems more appropriate

@tteschemacher is this decision documented somewhere?

In any case I think we should proceed with this PR since it also fixes some omp-issues (ping @loumalouomega)

[loumalouomega]

In a recent technical committee meeting we agreed to have the Is inside as a function including the projection. If the projection is done beforehand it was suggested to use IsInsideLocalSpace wouldn't it be the most consistent moving towards this interface.

Hm this sounds like a lot of work and is not applicable for e.g. solids
Now thinking about it the other approach seems more appropriate

@tteschemacher is this decision documented somewhere?

In any case I think we should proceed with this PR since it also fixes some omp-issues (ping @loumalouomega)

Why ping? I fixed already the issues?

[philbucher]

you are blocking

[loumalouomega]

Because I open #7734

[philbucher]

Because I open #7734

ah ok I missed that somehow

looks like too big/important to be merged before Christmas 😅

[loumalouomega]

Adding to @KratosMultiphysics/implementation-committee

[loumalouomega]

@philbucher ContactMechanicsApplication is failing despite the changes, I will need to check this, or even implement the method @KratosMultiphysics/implementation-committee suggested

[RiccardoRossi]

@loumalouomega @philbucher is this still active/relevant?

[loumalouomega]

Partially this is related wth #12637 (just line). We should agree how the IsInside should behave @KratosMultiphysics/technical-committee