The relationship between the GeometricField and the DimensionedField is a clear cut example of why inheriting for implementation does not work.
On one hand side, there is a public inheritance from DimensionedField:
template<class Type, template<class> class PatchField, class GeoMesh>
public DimensionedField<Type, GeoMesh>
and the OO principles state that the GeometricField IS-A DimensionedField, and should behave like that. As a result, following the Liskov substitution principle, every time a GeometricField is passed as a parameter to a function that takes a (const) reference of a DimensionedField, the call would work. It doesn't for assignment: the assignment operators have all been overriden
(overloaded for some types, but still they hide the assignment operators of the DimensionedField), in the GeometricField:
// Member operators
void operator=(const GeometricField<Type, PatchField, GeoMesh>&);
void operator=(const tmp<GeometricField<Type, PatchField, GeoMesh> >&);
void operator=(const dimensioned<Type>&);
void operator==(const tmp<GeometricField<Type, PatchField, GeoMesh> >&);
void operator==(const dimensioned<Type>&);
As a result, GeometricField is ... well... not working with DimensionedField. Such inheritance for implementation was maybe done on purpose, since the boundary fields are also dimensioned fields.
So if someone would code:
volScalarField p2 = p_1;
the semantical question would be "Should the internal field be assigned to, or one of the boundary fields?". To me, that question is retorical, the assignment operators should use those in DimensionedField, and the internal field should be assigned to. Besides that, 'operator==' was used precisely to introduce assignment to boundaries, so 'operator=' assignes to internal fields. The access to internal field is actually supported by the member function in GeometricField that gives a non-const access to the internal field:
//- Return dimensioned internal field
- breaks encapsulation
- forces the user to write the assignment by force
Also, you have placed the entire function implementation in a return statement. There is a reason functions have bodies. Put some code there. ;) Coding train wrecks like this one:
will make it very difficult for anyone trying to maintain your code to debug it and extend it.
Also, if you are working with temporary objects, take a look at how they are used in the code
, you can see how the interpolation/convection/diffusion schemes work, how they first initialize the temporary, operate on it, then return it. Make the code humanly readable, you will be grateful for it after enough time passes and you forget what you did - and so will all your users. :)
Here is a working example with that what is available in the GeometricField, off the top of my head. I am not sure if you can copy-paste-use that directly, but it is a start:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
const volScalarField& S2,
const volScalarField& F1
const Time& runTime = S2.time();
const fvMesh& mesh = S2.mesh();
// Initialize the temporary object that is to be returned.
tmp<volScalarField> resultTmp (
new volScalarField (
S2.dimensions(), // Judging by your max() call, you are returning S2 dimensions
// Get the non-const reference to the temporary object.
volScalarField& result = resultTmp();
// Keep dimensions of the result but re-name the field.
volScalarField h ("h", result);
// Get the non-const ref to the internal field which is of dimensioned field type.
DimensionedField<scalar, volMesh>& hInternal = h.dimensionedInternalField();
// Set the internal field values.
hInternal = pow(mesh.V(), 1.0 / 3.0);
// Calculate the result.
// FIXME: separate your huge return call into separate clearly readable calculations and put
// them here.
result = Foam::max(S2, h);
int main(int argc, char *argv)
tmp<volScalarField> result = Lvk(S2, F1);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nEnd\n" << endl;
Edit: also, try extracting your problem into a minimal working example, like the application I posted above, you'll understand your own problem better and get answers sooner. :)