CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   adding a C code to OpenFOAM (http://www.cfd-online.com/Forums/openfoam-programming-development/109549-adding-c-code-openfoam.html)

latvietis November 20, 2012 13:39

adding a C code to OpenFOAM
 
Greetings!

I have a code which I want to use in OpenFOAM. The idea is that at first there is given an initial value of "caurl" and then this loop goes trough all cells of my mesh and recalculates value of "caurl" which is afterwards written out and used again.

How it is possible to make this work with OpenFOAM? I have problems understanding all this class/object thing.

Lastly, b0, r0 and r2 are values taken from graph points and B2 is volScalarField.

The additional questions are:

1) Can I ignore the fact that B2 is defined as volScalarField and I'm trying to compare it with a float data type?

2) How it is possible to make this code work in every cell? Is it something similar to forAll(mesh.V(), celli)?

Code:

    float b0[n]={0,0.04,0.16,0.36,0.64,1.,1.44,1.69,1.96,2.25,2.56};
    float r0[n]={0.109E-03,0.109E-03,0.109E-03,0.11E-03,0.111E-03,0.112E-03,0.12E-03,0.143E-03,0.196E-03,0.417E-03,0.125E-02};
    float r2[n]={0,0,0,0,0.439E-05,0,0.248E-03,0.421E-03,0.466E-03,0.951E-02,0};

Code:

{
    float rval=1;
    for (int i=0; i<n; i++)
    {
        float* b0 = x;
        float* r0 = y;
        float* r2 = z;
        int kl,kr;

        if ( B2 <= b0[n-1])
            {
                kl = 0;
                kr = n-1;

                    if (kr - kl >= 1)
                    {
                        int k = (kr + kl)/2.;
                            if (b0[n-1] >= B2)
                                kr = k;
                            else
                                kl = k;
                    }
            float dx = b0[kr] - b0[kl];
            float du = (b0[kr] - B2)/dx;
            float dl = (B2 - b0[kl])/dx;
            float du2 = pow(du,2);
            float dl2 = pow(dl,2);
            rval = du * r0[kl] + dl * r0[kr] + ((du2 - 1.) * du * r2[kl] + (dl2 - 1.) * dl * r2[kr]) * pow(dx,2)/6.;
            }
            else
            {
            float dx = b0[n-1] - b0[n-2];
            float rder = -(r0[n-2] - r0[n-1])/dx + (+r2[n-2] + 2. * r2[n-1]) * dx/6.;
            rval = r0[n-1] + (B2 - b0[n-1]) * rder;
            }
    b0=NULL;
    r0=NULL;
    r2=NULL;
    }
    caurl=rval;
}

Sincerely,
Martin

akidess November 21, 2012 04:46

1) scalar and float are related, you only might get troubles with dimensions
2) Yes, it's exactly like that. forAll(B2, celli) will iterate through the volScalarField. If you use B2[celli], I think you will get a scalar value without dimensions.

- Anton


All times are GMT -4. The time now is 16:59.