### Converting Python code to C++ -- project function

226
views
0
13 months ago by

I have the following Python code

def ad(theta, a=0.1, Tr=5):
    return conditional(lt(theta,Tr), 0, a*(theta - Tr))

element = FiniteElement('CG', tetrahedron, 2)
V = FunctionSpace(mesh, element)
dt = 0.5
omega_n_1 = Function(V)
omega = Function(V)
theta = TrialFunction(V)
theta_n_1 = Function(V)

# cut code that solves a form for theta and assigns values to omega_n_1 and theta_n_1

omega.assign(project( (1 - (dt/2)*ad(theta_n_1))*omega_n_1 / (1 + (dt/2)*ad(theta)), V))

I tried to use this as equivalent C++ code

double ad(double theta_val, double a = 0.1, double Tr = 5)
{
    return theta_val < Tr ? 0.0 : a*(theta_val - Tr);
}

// and in a separate function
double dt = 0.5; auto theta = std::make_shared<dolfin::Function>(V); auto theta_n_1 = std::make_shared<dolfin::Function>(V);
auto omega = std::make_shared<dolfin::Function>(V);
auto omega_n_1 = std::make_shared<dolfin::Function>(V);// snipped code that puts values in the functions

for (std::size_t dof_num = 0; dof_num < omega->vector()->size(); ++dof_num)
{
    const double theta_n_1_val = theta_n_1->vector()->getitem(dof_num);
    const double theta_val = theta->vector()->getitem(dof_num);
    const double omega_n_1_val = omega_n_1->vector()->getitem(dof_num);
    const double omega_val =
        (1 - (dt / 2)*ad(theta_n_1_val))*omega_n_1_val / (1 + (dt / 2)*ad(theta_val));
    omega->vector()->setitem(dof_num, omega_val);
}

But these two do not return the same values for the same input

I am using fenics-2017-1

Thanks

Edited to highlight code snippets

Community: FEniCS Project