[CORE] Adding new RBF types in the RBFShapeFunctionsUtility class

[juancamarotti]

Description

This PR extends the RBFShapeFunctionsUtility class to include several additional Radial Basis Function (RBF) formulations commonly used for interpolation and mapping in multi-physics applications (e.g., FSI coupling).

Related PR
#13917

[matekelemen]

Idk for sure how often EvaluateRBF gets called, but I imagine it's very often, so it might make sense to make this efficient (and efficiently callable).

As it stands rn, there are (at least) two things that make this function more expensive than it needs to be:

1. It's defined in the source file. Yes, I usually bitch about too many definitions in the header, but this one really belongs there to give compilers a chance to inline it.

2. Now there's a switch statement in the function that

   • makes the function larger, so it's less likely to be inlined
   • adds (a possibly expensive) branching

Instead of passing an enum that represents an operation every single time that operation is carried out, I'd rather pass it outside the loop, and only pass the (inlinable) operation itself to this function.

Based on our earlier chats, I'm not going to try to further explain myself. I'll just show an analogy of the current code and what I'd want to see.

Current code

enum class OpType {
    Op1,
    Op2
};

double op1(double in) {
    double out;
    // compute some stuff ...
    return out;
}

double op2(double in) {
    double out;
    // compute some other stuff
    return out;
}

void applyOp(const double* begin,
             const double* end,
             double* out,
             OpType opType)
{
    for (auto in=begin; in!=end; ++in, ++out) {
        switch (opType) {
            case OpType::Op1: {
                // Do some stuff on the component.
                *out = op1(*in);
                break;
            }
            case OpType::Op2: {
                // Do some other stuff on component.
                *out = op2(*in);
                break;
            }
            //default: throw an exception
        }
    }
}

int main() {
    double in[100], out[100];
    
    OpType opType = OpType::Op2;

    applyOp(in, in + 100, out, opType);
}

What I want to see

enum class OpType {
    Op1,
    Op2
};

double op1(double in) {
    double out;
    // compute some stuff ...
    return out;
}

double op2(double in) {
    double out;
    // compute some other stuff
    return out;
}

template <class TOp>
void applyOp(const double* begin,
             const double* end,
             double* out,
             TOp&& op)
{
    for (auto in=begin; in!=end; ++in, ++out) 
        *out = op(*in); 
}

int main() {
    double in[100], out[100];
    
    OpType op_type = OpType::Op2;

    switch (op_type) {
        case OpType::Op1: {
            // Do some stuff on the array.
            applyOp(in, in + 100, out, op1);
            break;
        }
        case OpType::Op2: {
            // Do some other stuff on the array.
            applyOp(in, in + 100, out, op2);
            break;
        }
        //default: throw an exception
    }
}

[juancamarotti]

Ok, I tried to do this in the code, with the only difference that I did not implement theapplyOp()templated function as I think it is not necessary in my code (because I never call this function on arrays of doubles but on doubles).

[juancamarotti]

Then, inside the Mapper, I can have a switch() statement to decide which function to call.

[rubenzorrilla]

Which is the point on adding new RBF functionals if then you hardcode current behavior?

[juancamarotti]

We need these new functions for the RBF mapper, you can check my previous PR regarding this...

[rubenzorrilla]

Fine, but how are you using them? If these are not used through this class API which is the point of adding them in here? My point is that you can make them optional also here rather than hardcoding current behavior.

[juancamarotti]

I wanted to make my changes as non-intrusive as possible. Therefore, I hardcoded the current behavior to avoid modifying all the files where it's used. Would you like me to proceed the way you commented before?

[juancamarotti]

What do you think @matekelemen about the new implementation?

[philbucher]

LGTM

[juancamarotti]

@rubenzorrilla do you agree with these changes now?

[rubenzorrilla]

@rubenzorrilla do you agree with these changes now?

No as you're breaking current API besides making it not consistent with the MLSShapeFunctions utility, something that was done at purpose. Use default arguments, as we do w/ the QR in order to avoid so. As an alternative you can also overload the method.

[rubenzorrilla]

Blocking as it breaks current API besides making it not consistent with the homonymous MLS utility.

[juancamarotti]

@rubenzorrilla do you agree with these changes now?

No as you're breaking current API besides making it not consistent with the MLSShapeFunctions utility, something that was done at purpose. Use default arguments, as we do w/ the QR in order to avoid so. As an alternative you can also overload the method.

Could you elaborate a little bit more on this comment? I don't quite get your point regarding breaking the current API...

[juancamarotti]

@philbucher @matekelemen What do you think about adding a utility within the MappingApplication for RBF evaluation? We're using it in a way that's quite different from the core implementation, and in my opinion, it wouldn’t make sense to change the existing core API for this purpose.

[juancamarotti]

I will close this PR and move my implementation to an utility inside the MappingApplication

FYI @philbucher