How to Implement a “Function Object” suitable for NDSolve and Finite Element Method?

I have a bunch of methods that shoud return an expression that behave like a function, i.e. an expression that, when supplied with suitable arguments, evaluate to something.

A simple Function is not able to evaluate only for certain type of arguments, for example only for numeric arguments, or evaluate differently based on the number or type of arguments.

A possible way is to return an Module symbol with some DownValues attached, but, I don't like so much to return something like name\$nnnn and moreover, from what I understand, these kind of symbols should be explicitely ClearAll once no more needed to avoid memory leaks.

Instead of defining a function object for each method or each kind of return type, I tried to define a common FunctionObject that I can return from every method but I encountered some problem whan passing these objects as a source function to NDSolve.

The obvous way with DownValues. Everything works as expected.

Clear[tfn]
tfn[p__?NumericQ] := {1, 1, 1} Norm[{p}]

tfn[x, y, z]
tfn[1, 2, 3]

NDSolveValue[{Laplacian[u[x, y, z], {x, y, z}] ==
Indexed[tfn[x, y, z], 1],
DirichletCondition[u[x, y, z] == 0, True]}, u, {x, y, z} \[Element] Ball[]]


tfn[x, y, z]

{Sqrt[14], Sqrt[14], Sqrt[14]}

InterpolatingFunction[...]

The first tentative with a "wrapper" for a pure Function. Everything apparently works as expected.

Clear[tfn]
FunctionPattern[fn_Function][args___] /;
MatchQ[HoldComplete[args], HoldComplete[___?NumericQ]] := fn@args
tfn = FunctionPattern[{1, 1, 1} Norm@{##} &];

tfn[x, y, z]
tfn[1, 2, 3]

NDSolveValue[{Laplacian[u[x, y, z], {x, y, z}] ==
Indexed[tfn[x, y, z], 1],
DirichletCondition[u[x, y, z] == 0, True]}, u, {x, y, z} \[Element] Ball[]]


FunctionPattern[{1, 1, 1} Norm[{##1}] &][x, y, z]

{Sqrt[14], Sqrt[14], Sqrt[14]}

InterpolatingFunction[{{-0.999955, 0.999955}, {-1., 1.}, {-1., 1.}}, <>]

A more general way. As stated, I want also to be able to handle differently different types of arguments in a more general way. I tried the following approach. Really, before the following, I tried alternative and possibly more efficient definitions, but the following is the more similar to the previous working definition:

ClearAll[FunctionObject2]
FunctionObject2[defs : {(_Rule | _RuleDelayed) ...}][args___] /;
AnyTrue[defs,
With[{form = First@#},
MatchQ[HoldComplete[{args}], HoldComplete[form]]] &] :=
Replace[{args}, defs]

Clear[tfn]
tfn = FunctionObject2[{{p__?NumericQ} :> {1, 1, 1} Norm[{p}]}];

tfn[x, y, z]
tfn[1, 2, 3]


FunctionObject2[{{p__?NumericQ} :> {1, 1, 1} Norm[{p}]}][x, y, z]

{Sqrt[14], Sqrt[14], Sqrt[14]}

This is good, but when suppling this function object to NDSolve, NDSolve complains (but still provide a result in this simplified scenario):

NDSolveValue[{Laplacian[u[x, y, z], {x, y, z}] ==
Indexed[tfn[x, y, z], 1],
DirichletCondition[u[x, y, z] == 0, True]}, u, {x, y, z} \[Element] Ball[]]


During evaluation of In[426]:= CompiledFunction::cfex: Could not complete external evaluation at instruction 11; proceeding with uncompiled evaluation. >>

During evaluation of In[426]:= NDSolveValue::femdpop: The FEMLoadElements operator failed. >>

During evaluation of In[426]:= NDSolveValue::femdpop: The FEMLoadElements operator failed. >>

InterpolatingFunction[{{-0.999955, 0.999955}, {-1., 1.}, {-1., 1.}}, <>]

In my real, more complex, scenario I get different messages and no result is provided:

NDSolveValue::femcnsd: The PDE coefficient [...] FunctionObject[...][3.6358*10^9,6.37814*10^6 x, 6.37814*10^6 y, z] does not evaluate to a numeric scalar. >>

What can be the deifference, from the point of view of evaluation, between the second and third way, that can cuase the third to fail?

Any help is appreciated. Any better design too.