Skip to content

GitLab

Commit bee953fd authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

remove functionPython and functionNumpy

parent bd9a0593
Hide whitespace changes
Inline Side-by-side
dgpy/scripts/levelset.py
View file @ bee953fd
......@@ -7,13 +7,9 @@ import os
#-------------------------------------------------
#-- Generalfunctions
#-------------------------------------------------
def areaF(FCT, sol):
for i in range(0,FCT.size1()):
phi = sol.get(i,0)
value = 0.0
if (phi < 0.0):
value = 1.0
FCT.set(i,0, value)
def areaF(cm, sol):
phi = cm.get(sol)[:,[0]]
return numpy.where(phi<0, 1., 0)
#-------------------------------------------------
#-- CLASS DEFINITION
......@@ -154,7 +150,7 @@ class levelset(object):
def computeArea(self, t):
fun = open("Area.dat","a")
self.AREA = functionPython(1, areaF, [self.solutions[0].getFunction()])
self.AREA = functorNumpy(lambda cm : areaF(cm, self.solutions[0].getFunction()))
integ=dgFunctionIntegrator(self.groups[0], self.AREA)
localArea=fullMatrixDouble(1,1)
integ.compute(localArea)
......
dgpy/scripts/slim3d_private.py
View file @ bee953fd
......@@ -5,6 +5,7 @@ import dgpy
from dgpy.scripts import slim_private
import numpy
def slim3d_setup(loop):
slimSolver = loop._slimSolver
eq = loop._equations
......@@ -31,21 +32,20 @@ def slim3d_setup(loop):
f = slimSolver.functions
f.initializeBath()
def mergeUVNumpy(cmap, val, u, v) :
val[:,0] = u[:]
val[:,1] = v[:]
def mergeUVNumpy(cm, u, v) :
return numpy.hstack([cm.get(u)[:,[0]], cm.get(v)[:,[0]]])
if eq._wind_stress:
if eq._wind_stress[0] == 'speed':
eq._wind_speed_x = slim_private._load_function(eq._wind_stress[1], slimSolver.groups2d)
eq._wind_speed_y = slim_private._load_function(eq._wind_stress[2], slimSolver.groups2d)
f.windFunc = dgpy.functionNumpy(2, mergeUVNumpy, [eq._wind_speed_x, eq._wind_speed_y])
f.windFunc = functorNumpy(lambda cm : mergeUVNumpy(cm, eq.wind_speed_x, eq.wind_speed_y))
f.windStressFunc = dgpy.windCubicStress(f.windFunc, eq._wind_stress[3])
# TODO Fatal if no wind speed and solving sediment
elif eq._wind_stress[0] == 'stress':
eq._wind_stress_x = slim_private._load_function(eq._wind_stress[1], slimSolver.groups2d)
eq._wind_stress_y = slim_private._load_function(eq._wind_stress[2], slimSolver.groups2d)
f.windStressFunc = dgpy.functionNumpy(2, mergeUVNumpy, [eq._wind_stress_x, eq._wind_stress_y])
f.windStressFunc = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm,eq._wind_stress_x, eq._wind_stress_y))
if eq._initial_elevation:
f.etaInitFunc = slim_private._load_function(eq._initial_elevation, slimSolver.groups2d)
......@@ -56,9 +56,7 @@ def slim3d_setup(loop):
if eq._initial_salinity[0] == 'netcdf':
f.SInitFunc = slim_private._load_function(eq._initial_salinity[1], slimSolver.groups3d)
elif eq._initial_salinity[0] == 'vertical_gradient':
def zFunc(cmap, val, xyz) :
val[:] = eq._initial_salinity[2] + xyz[:,2] * eq._initial_salinity[3]
f.SInitFunc = dgpy.functionNumpy(1, zFunc, [d.xyzOrigDof.getFunction()])
f.SInitFunc = dgpy.functorNumpy(lambda cm : eq._initial_salinity[2]+cm.get(d.xyzOrigDof.getFunction())[:,[2]]*eq._initial_salinity[3])
else:
dgpy.Msg.Fatal('Unknown mode for initial salinity:' + eq._initial_salinity[0])
elif (not eq._linear_density) or (eq._linear_density[0] == 'salinity'):
......@@ -68,9 +66,7 @@ def slim3d_setup(loop):
if eq._initial_temperature[0] == 'netcdf':
f.TInitFunc = slim_private._load_function(eq._initial_temperature[1], slimSolver.groups3d)
elif eq._initial_temperature[0] == 'vertical_gradient':
def zFunc(cmap, val, xyz) :
val[:] = eq._initial_temperature[2] + xyz[:,2] * eq._initial_temperature[3]
f.TInitFunc = dgpy.functionNumpy(1, zFunc, [d.xyzOrigDof.getFunction()])
f.TInitFunc = dgpy.functorNumpy(lambda cm : eq._initial_temperature[2]+cm.get(d.xyzOrigDof.getFunction())[:,[2]]*eq._initial_temperature[3])
else:
dgpy.Msg.Fatal('Unknown mode for initial temperature:' + eq._initial_temperature[0])
elif (not eq._linear_density) or (eq._linear_density[0] == 'temperature'):
......@@ -157,23 +153,22 @@ def slim3d_setup(loop):
g = eq._gravity
# Open Boundaries
def uvOpenGivenEtaNumpy(cmap, val, uv, normals) :
un = uv[:,0] * normals[:,0] + uv[:,1] * normals[:,1]
val[:,0] = un * normals[:,0]
val[:,1] = un * normals[:,1]
def uvOpenFreeNumpy(cmap, val, bath, eta, normals) :
val[:,0] = numpy.sqrt(g / (bath[:] + eta[:]) ) * eta[:] * normals[:,0]
val[:,1] = numpy.sqrt(g / (bath[:] + eta[:]) ) * eta[:] * normals[:,1]
def etaOpenGivenUVNumpy(cmap, val, uv, bath, eta, normals) :
un = uv[:,0] * normals[:,0] + uv[:,1] * normals[:,1]
val[:] = numpy.sqrt( (bath[:]+eta[:]) / g) * un
def transport2UVNumpy(cmap, val, uv, bath, eta) :
val[:,0] = uv[:,0] / ( bath[:] + eta[:] )
val[:,1] = uv[:,1] / ( bath[:] + eta[:] )
def mergeEtaUV(cmap, val, eta, uv):
val[:,0] = eta[:]
val[:,1] = uv[:,0]
val[:,2] = uv[:,1]
def uvOpenGivenEtaNumpy(cm, u, normals) :
u = cm.get(uv)
n = cm.get(normals)
un = u[:,[0]] * n[:,[0]] + u[:,[1]] * n[:,[1]]
return un * n
def uvOpenFreeNumpy(cm, bath, eta, normals) :
c = numpy.sqrt(g/(cm.get(bath)+cm.get(eta)))
return c*cm.get(eta)*cm.get(normals)
def etaOpenGivenUVNumpy(cm, u, bath, eta, normals) :
u = cm.get(uv)
n = cm.get(normals)
un = u[:,[0]] * n[:,[0]] + u[:,[1]] * n[:,[1]]
return numpy.sqrt( (cm.get(bath)+cm.get(eta)) / g) * un
def transport2UVNumpy(cm, uv, bath, eta) :
return cm.get(uv) / (cm.get(bath) + cm.get(eta))
eta = f.eta2dFunc
uv = d.uvDof.getFunction()
......@@ -203,13 +198,13 @@ def slim3d_setup(loop):
eq._uAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d)
eq._vAv2d_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups2d)
if openBnd.transport:
eq._uvTransport_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uv_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvTransport_open[openBnd], f.bathFunc2d, eta])
eq._uvAv2dTransport_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvAv2dTransport_open[openBnd], f.bathFunc2d, eta])
eq._uvTransport_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._u_open[openBnd], eq._v_open[openBnd]))
eq._uv_open[openBnd] = dgpy.functorNumpy(lambda cm : transport2UVNumpy(cm, eq._uvTransport_open[openBnd], f.bathFunc2d, eta))
eq._uvAv2dTransport_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]))
eq._uvAv2d_open[openBnd] = dgpy.functorNumpy(lambda cm : transport2UVNumpy(cm, eq._uvAv2dTransport_open[openBnd], f.bathFunc2d, eta))
else:
eq._uv_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._uv_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._u_open[openBnd], eq._v_open[openBnd]))
eq._uvAv2d_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]))
eq._eta_open[openBnd] = slim_private._load_function(openBnd.eta, slimSolver.groups2d)
elif openBnd.u:
eq._u_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d)
......@@ -217,17 +212,17 @@ def slim3d_setup(loop):
eq._uAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d)
eq._vAv2d_open[openBnd] = slim_private._load_function(openBnd.v, slimSolver.groups2d)
if openBnd.transport:
eq._uvTransport_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uv_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvTransport_open[openBnd], f.bathFunc2d, eta])
eq._uvAv2dTransport_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvAv2dTransport_open[openBnd], f.bathFunc2d, eta])
eq._uvTransport_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._u_open[openBnd], eq._v_open[openBnd]))
eq._uv_open[openBnd] = dgpy.functorNumpy(lambda cm : transport2UVNumpy(cm, eq._uvTransport_open[openBnd], f.bathFunc2d, eta))
eq._uvAv2dTransport_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]))
eq._uvAv2d_open[openBnd] = dgpy.functorNumpy(lambda cm : transport2UVNumpy(cm, eq._uvAv2dTransport_open[openBnd], f.bathFunc2d, eta))
else:
eq._uv_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._eta_open[openBnd] = dgpy.functionNumpy(1, etaOpenGivenUVNumpy, [uvAv2d, f.bathFunc2d, eta, dgpy.function.getNormals()])
eq._uv_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._u_open[openBnd], eq._v_open[openBnd]))
eq._uvAv2d_open[openBnd] = dgpy.functorNumpy(lambda cm : mergeUVNumpy(cm, eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]))
eq._eta_open[openBnd] = dgpy.functorNumpy(lambda cm : etaOpenGivenUVNumpy(cm, uvAv2d, f.bathFunc2d, eta, dgpy.function.getNormals()))
elif openBnd.eta:
eq._uv_open[openBnd] = dgpy.functionNumpy(2, uvOpenGivenEtaNumpy, [uv, dgpy.function.getNormals()])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, uvOpenGivenEtaNumpy, [uvAv2d, dgpy.function.getNormals()])
eq._uv_open[openBnd] = dgpy.functorNumpy(lambda cm : uvOpenGivenEtaNumpy(cm, uv, dgpy.function.getNormals()))
eq._uvAv2d_open[openBnd] = dgpy.functorNumpy(lambda cm : uvOpenGivenEtaNumpy(cm, uvAv2d, dgpy.function.getNormals()))
eq._eta_open[openBnd] = slim_private._load_function(openBnd.eta, slimSolver.groups2d)
elif openBnd.flux:
integInterface = dgpy.dgFunctionIntegratorInterface(slimSolver.groups2d, f.bathFunc2d)
......@@ -240,8 +235,8 @@ def slim3d_setup(loop):
eq._uvAv2d_open[openBnd] = dgpy.slimFlowRateToVelocity(slimSolver.groups2d, eq._flux_section[openBnd], eq._flux2d[openBnd], f.bathFunc2d, eta)
eq._eta_open[openBnd] = eta
else:
eq._uv_open[openBnd] = dgpy.functionNumpy(2, uvOpenFreeNumpy, [f.bathFunc2d, eta, dgpy.function.getNormals()])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, uvOpenFreeNumpy, [f.bathFunc2d, eta, dgpy.function.getNormals()])
eq._uv_open[openBnd] = dgpy.functorNumpy(lambda cm : uvOpenFreeNumpy(cm, f.bathFunc2d, eta, dgpy.function.getNormals()))
eq._uvAv2d_open[openBnd] = dgpy.functorNumpy(lambda cm : uvOpenFreeNumpy(cm, f.bathFunc2d, eta, dgpy.function.getNormals()))
eq._eta_open[openBnd] = eta
if openBnd.salinity:
if openBnd.salinity == 'vertical_gradient':
......@@ -253,7 +248,7 @@ def slim3d_setup(loop):
eq._T_open[openBnd] = f.TInitFunc
else:
eq._T_open[openBnd] = slim_private._load_function(openBnd.temperature, slimSolver.groups3d)
eq._sw2D_open[openBnd] = dgpy.functionNumpy(3, mergeEtaUV, [eq._eta_open[openBnd], eq._uvAv2d_open[openBnd]])
eq._sw2D_open[openBnd] = dgpy.functorNumpy(lambda cm : numpy.hstack([cm.get(eq._eta_open[openBnd]), cm.get(eq._uvAv2d_open[openBnd])]))
sw2DEq = e.sw2DEq
sw2DEq.setFrom3D(True)
......
dgpy/scripts/slimPre.py
View file @ bee953fd
......@@ -577,9 +577,11 @@ def residual_flow(region, flow, bathymetry, data_file_name=None):
ny = ny/norm
bathFunc = slim_private._load_function(bathymetry, groups)
def orientedSurfNumpy(cmap, val, h, n) :
val[:] = h[:] * (n[:,0]*nx + n[:,1]*ny)
f = dgpy.functionNumpy(1, orientedSurfNumpy, [bathFunc, dgpy.function.getNormals()])
def orientedSurfNumpy(cm) :
n = cm.get(dgpy.function.getNormals())
h = cm.get(bathFunc)
return h[:,[0]] * (n[:,[0]]*nx + n[:,[1]]*ny)
f = dgpy.functorNumpy(orientedSurfNumpy)
integInterface = dgpy.dgFunctionIntegratorInterface(groups, f)
section = 0
......@@ -589,10 +591,10 @@ def residual_flow(region, flow, bathymetry, data_file_name=None):
integInterface.compute(tag, sectionFM)
section += sectionFM.get(0,0)
def transportNumpy(cmap, val, h) : # eta is neglected here
val[:,0] = -h[:] * flow / section * nx
val[:,1] = -h[:] * flow / section * ny
f = dgpy.functionNumpy(2, transportNumpy, [bathFunc])
def transportNumpy(cm) : # eta is neglected here
f = cm.get(bathFunc)*flow/section
return numpy.hstack([-f*nx,-f*ny])
f = dgpy.functorNumpy(transportNumpy)
vals = region._evaluateFunctor(f, 2)
if region._empty:
vals = np.empty((1,2))
......
dgpy/scripts/spectralAnalysis1D.py
View file @ bee953fd
......@@ -173,11 +173,8 @@ def spectralAnalysis1D(LAW, FLUX, ORDER, N, NFFT,
print("--- Saving eigenvectors and computing wavenumbers")
eigenVecRealDof = dgDofContainer(groups,nFields)
eigenVecImagDof = dgDofContainer(groups,nFields)
def eigenVecFuncTMP(val,func):
for i in range(0,func.size1()):
val.set(i,0,func.get(i,0)) # For the first field only
eigenVecRealFunc = functionPython(1,eigenVecFuncTMP,[eigenVecRealDof.getFunction()])
eigenVecImagFunc = functionPython(1,eigenVecFuncTMP,[eigenVecImagDof.getFunction()])
eigenVecRealFunc = functorNumpy(lambda cm : cm.get(eigenVecRealDof.getFunction())[:,[0]]) # For the first field only
eigenVecImagFunc = functorNumpy(lambda cm : cm.get(eigenVecImagDof.getFunction())[:,[0]]) # For the first field only
eigenVecRealFuncEvaluator = dgFunctionEvaluator(groups, eigenVecRealFunc)
eigenVecImagFuncEvaluator = dgFunctionEvaluator(groups, eigenVecImagFunc)
valueImag = fullMatrixDouble(1,1)
......
dgpy/scripts/spectralAnalysis2D.py
View file @ bee953fd
......@@ -214,11 +214,8 @@ def spectralAnalysis2D(LAW, FLUX, THETA, MESH, ORDER, N, NFFT,
print("--- Saving eigenvectors and computing wavenumbers")
eigenVecRealDof = dgDofContainer(groups,nFields)
eigenVecImagDof = dgDofContainer(groups,nFields)
def eigenVecFuncTMP(val,func):
for i in range(0,func.size1()):
val.set(i,0,func.get(i,0)) # For the first field only
eigenVecRealFunc = functionPython(1,eigenVecFuncTMP,[eigenVecRealDof.getFunction()])
eigenVecImagFunc = functionPython(1,eigenVecFuncTMP,[eigenVecImagDof.getFunction()])
eigenVecRealFunc = functorNumpy(lambda cm : cm.get(eigenVecRealDof.getFunction())[:,[0]]) # For the first field only
eigenVecImagFunc = functorNumpy(lambda cm : cm.get(eigenVecImagDof.getFunction())[:,[0]]) # For the first field only
eigenVecRealFuncEvaluator = dgFunctionEvaluator(groups, eigenVecRealFunc)
eigenVecImagFuncEvaluator = dgFunctionEvaluator(groups, eigenVecImagFunc)
valueReal = fullMatrixDouble(1,1)
......
function/dgFunction.i
View file @ bee953fd
......@@ -10,9 +10,7 @@
#include "functor.h"
#include "function.h"
#include "functionGeneric.h"
#include "functionPython.h"
#include "functionNumpy.h"
#include "dgFunctionIntegrator.h"
#include "dgFunctionIntegratorInterface.h"
#include "dgFunctionEvaluator.h"
......@@ -42,7 +40,6 @@ namespace std {
%include "function.h"
%warnfilter(509) functionConstantByTag;
%include "functionGeneric.h"
%include "functionPython.h"
%include "functionNumpy.h"
%include "dgFunctionIntegrator.h"
%include "dgFunctionIntegratorInterface.h"
......
function/functionNumpy.h
View file @ bee953fd
......@@ -10,76 +10,6 @@
#endif
#include "functor.h"
class functionNumpy : public functor {
PyObject *_pycallback;
int _nbCol;
std::vector<std::pair<const functor*, int> > _deps;
static PyObject *M2A(const fullMatrix<double> *m) {
npy_intp n[2] = {m->size1(), m->size2()};
return PyArray_New(&PyArray_Type, m->size2() > 1 ? 2 : 1, n, NPY_DOUBLE, NULL, (void*)m->getDataPtr(), 0, NPY_ARRAY_FARRAY, NULL);
}
static void _checkImported() {
static bool _arrayImported = false;
if (! _arrayImported){
_import_array();
_arrayImported = true;
}
}
public:
void operator()(functorCache &m, fullMatrix<double> &res)const
{
PyObject *pyargs;
PyObject *cmpy = SWIG_NewPointerObj((void*) &m, SWIGTYPE_p_dataCacheMap, 0);
res.resize(m.nEvaluationPoint(), _nbCol);
std::vector<const fullMatrix<double>*> args;
for (auto &d : _deps)
args.push_back(&m.get(*d.first, d.second));
switch(args.size()) {
case 0 : pyargs = Py_BuildValue("(NN)", cmpy, M2A(&res)); break;
case 1 : pyargs = Py_BuildValue("(NNN)", cmpy, M2A(&res), M2A(args[0])); break;
case 2 : pyargs = Py_BuildValue("(NNNN)", cmpy, M2A(&res), M2A(args[0]), M2A(args[1])); break;
case 3 : pyargs = Py_BuildValue("(NNNNN)", cmpy, M2A(&res), M2A(args[0]), M2A(args[1]), M2A(args[2])); break;
case 4 : pyargs = Py_BuildValue("(NNNNNN)", cmpy, M2A(&res), M2A(args[0]), M2A(args[1]), M2A(args[2]), M2A(args[3])); break;
case 5 : pyargs = Py_BuildValue("(NNNNNNN)", cmpy, M2A(&res), M2A(args[0]), M2A(args[1]), M2A(args[2]), M2A(args[3]), M2A(args[4])); break;
case 6 : pyargs = Py_BuildValue("(NNNNNNNN)", cmpy, M2A(&res), M2A(args[0]), M2A(args[1]), M2A(args[2]), M2A(args[3]), M2A(args[4]), M2A(args[5])); break;
case 7 : pyargs = Py_BuildValue("(NNNNNNNNN)", cmpy, M2A(&res), M2A(args[0]), M2A(args[1]), M2A(args[2]), M2A(args[3]), M2A(args[4]), M2A(args[5]), M2A(args[6])); break;
case 8 : pyargs = Py_BuildValue("(NNNNNNNNNN)", cmpy, M2A(&res), M2A(args[0]), M2A(args[1]), M2A(args[2]), M2A(args[3]), M2A(args[4]), M2A(args[5]), M2A(args[6]), M2A(args[7])); break;
default:Msg::Error("python function not implemented for more than 8 arguments");
}
PyObject *result = PyEval_CallObject(_pycallback, pyargs);
if (result) {
Py_DECREF(result);
}
else {
PyErr_Print();
Msg::Fatal("An error occurs in the python function.");
}
Py_DECREF(pyargs);
}
functionNumpy (int nbCol, PyObject *callback, std::vector<const functor*> dependencies)
: _pycallback(callback)
{
_nbCol = nbCol;
Py_INCREF(_pycallback);
_deps.reserve(dependencies.size());
for (auto &d : dependencies) {
_deps.push_back(std::make_pair(d,-1));
}
_checkImported();
}
functionNumpy (int nbCol, PyObject *callback, std::vector<std::pair<const functor*, int> > dependencies)
: _pycallback(callback)
{
_nbCol = nbCol;
Py_INCREF(_pycallback);
_deps = dependencies;
_checkImported();
}
~functionNumpy() {
Py_DECREF(_pycallback);
}
};
class functorNumpy : public functor {
PyObject *_pycallback;
......
function/tests/functionNumpy/test.py
View file @ bee953fd
......@@ -11,29 +11,18 @@ wind = {}
XYZ = function.getCoordinates()
#### NUMPY
def coriolisNumpy(cmap, value, xyz):
value[:] = 1e-4 + 2e-11 * xyz[:,1]
def coriolisNumpy(cm) :
y = cm.get(XYZ)[:,[1]]
return 1e-4 + 2e-11 * y
def windNumpy(cmap, value, xyz):
value[:,0] = numpy.sin(xyz[:,1]/1e6)/1e6
value[:,1] = 0.
def windNumpy(cm):
v = numpy.zeros((cm.nEvaluationPoint(), 2))
y = cm.get(XYZ)[:,1]
v[:,0] = numpy.sin(y/1e6)/1e6
return v
coriolis["NUMPY"] = functionNumpy(1, coriolisNumpy, [XYZ])
wind["NUMPY"] = functionNumpy(2, windNumpy, [XYZ])
#### PYTHON
def coriolisPython(value, xyz):
for i in range(xyz.size1()):
value.set(i, 0, 1e-4 + 2e-11 * xyz(i,1))
def windPython(value, xyz):
for i in range(xyz.size1()):
value.set(i, 0, math.sin(xyz(i,1)/1e6)/1e6)
value.set(i, 1, 0.)
coriolis["PYTHON"] = functionPython(1, coriolisPython, [XYZ])
wind["PYTHON"] = functionPython(2, windPython, [XYZ])
coriolis["NUMPY"] = functorNumpy(coriolisNumpy)
wind["NUMPY"] = functorNumpy(windNumpy)
### C
CCode = """
......@@ -67,7 +56,7 @@ groups = dgGroupCollection("stommel_square.msh")
order = 1
results = {}
print("%10s %10s %20s" % ("Mode", "Time", "Solution norm"))
for mode in ["C", "NUMPY", "PYTHON"] :
for mode in ["C", "NUMPY"] :
claw = dgConservationLawShallowWater2d()
solution = dgDofContainer(groups, claw.getNbFields())
solution.setAll(0.)
......
linearSystem/linearSystem.h
View file @ bee953fd
......@@ -25,7 +25,7 @@ class linearSystemBase {
virtual void zeroSolution() = 0;
virtual int systemSolve() = 0;
// x = A*b
virtual int matMult() { return 0; }
virtual int matMult() {return 0; }
void setParameter (std::string key, std::string value);
std::string getParameter(std::string key) const;
......
linearSystem/linearSystemPETSc.h
View file @ bee953fd
......@@ -148,7 +148,7 @@ class linearSystemPETSc : public linearSystem<scalar> {
void printMatlab(const char *filename) const{};
virtual int systemSolve() { return 0; }
double normInfRightHandSide() const{return 0;}
virtual int matMult() { return 0; }
virtual int matMult() {return 0; }
};
#endif
#endif
mesh/mesh.i
View file @ bee953fd
......@@ -114,7 +114,7 @@
}
return ll;
}
fullMatrix<double> stereo2CartMat(const fullMatrix<double> &stereo, dataCacheMap *cacheMap){
fullMatrix<double> stereo2CartMat(const fullMatrix<double> &stereo, functorCache *cacheMap){
fullMatrix<double> xyz(stereo.size1(),stereo.size2());
for (int i=0; i<stereo.size1();i++){
$self->stereo2Cart(stereo(i,0), stereo(i,1), stereo(i,2), cacheMap, xyz(i,0), xyz(i,1), xyz(i,2));
......@@ -128,14 +128,14 @@
}
return vc;
}
fullMatrix<double> cartToStereoVecMat(const fullMatrix<double> &stereo, const fullMatrix<double> &v, dataCacheMap *cacheMap){
fullMatrix<double> cartToStereoVecMat(const fullMatrix<double> &stereo, const fullMatrix<double> &v, functorCache *cacheMap){
fullMatrix<double> vs(v.size1(),v.size2());
for (int i=0; i<stereo.size1();i++){
$self->cartToStereoVec(stereo(i,0), stereo(i,1), stereo(i,2), v(i,0), v(i,1), v(i,2), cacheMap, vs(i,0), vs(i,1), vs(i,2));
}
return vs;
}
fullMatrix<double> stereoToCartVecMat(const fullMatrix<double> &stereo, const fullMatrix<double> &v, dataCacheMap *cacheMap){
fullMatrix<double> stereoToCartVecMat(const fullMatrix<double> &stereo, const fullMatrix<double> &v, functorCache *cacheMap){
fullMatrix<double> vc(v.size1(),v.size2());
for (int i=0; i<stereo.size1();i++){
$self->stereoToCartVec(stereo(i,0), stereo(i,1), stereo(i,2), v(i,0), v(i,1), v(i,2), cacheMap, vc(i,0), vc(i,1), vc(i,2));
......
mesh/tests/CircleTransform/test.py
View file @ bee953fd
from dgpy.scripts.extrude import *
from dgpy import *
from math import *
import numpy
import time, os, sys
from partitionThenExtrude import *
from dgpy.scripts import Common
......@@ -35,24 +36,23 @@ for i in range(4):
XYZ = groups.getFunctionCoordinates()
def toIntegrate(FCT, XYZ) :
for i in range (0,XYZ.size1()) :
FCT.set(i,0,1)
def toIntegrate(cm) :
return numpy.ones((cm.nEvaluationPoint(), 1))
def toIntegrate2(FCT, XYZ) :
for i in range (0,XYZ.size1()) :
FCT.set(i,0,sqrt(XYZ(i,0)**2+XYZ(i,1)**2))
def toIntegrate2(cm) :
xyz = cm.get(xyzCircleInteg)
return numpy.hypot(xyz[:,[0]], xyz[:,[1]])
xyzCircle = transformToCircle(circleTransform)
xyzCircleInteg = functionPrecomputed(groups, integOrder)
xyzCircleInteg.compute(xyzCircle)
toInteg=functionPython(1, toIntegrate, [xyzCircleInteg])
toInteg=functorNumpy(toIntegrate)
integrator = dgFunctionIntegratorInterface(groups, toInteg)
integral = fullMatrixDouble(1,1)
integrator.compute("Seam",integral)
toInteg2=functionPython(1, toIntegrate2, [xyzCircleInteg])
toInteg2=functorNumpy(toIntegrate2)
integrator2 = dgFunctionIntegrator(groups, toInteg2)
integral2 = fullMatrixDouble(1,1)
integrator2.compute(integral2)
......@@ -72,18 +72,12 @@ for i in range(4):
integrals[i].append(integral4(0,0))
def exportFct(FCT,XYZ,xyzCircle):
for i in range(0,FCT.size1()):
FCT.set(i,0,XYZ(i,0))
FCT.set(i,1,XYZ(i,1))
FCT.set(i,2,XYZ(i,2))
FCT.set(i,3,xyzCircle(i,0))
FCT.set(i,4,xyzCircle(i,1))
FCT.set(i,5,xyzCircle(i,2))
def exportFct(cm) :
return cm.get(XYZ) + cm.get(xyzCircle)
nCompExp=[3,3]
namesExp=["xyz","xyzCircle"]
exportFunction=functionPython(sum(nCompExp), exportFct, [XYZ, xyzCircle])
exportFunction=functorNumpy(exportFct)
dof.exportFunctionVtk(exportFunction,'output', 0, 0,"solution",nCompExp,namesExp, xyzCircle)
print ("-------------------------------------------------------")
......
mesh/tests/periodicAdvection/periodicAdvection.py
View file @ bee953fd
# Import for python
from dgpy import *
from math import *
import numpy
import time, os
......@@ -24,25 +25,23 @@ Vz=0. # Advective Speed in z-direction
Function for Initial Condition
"""
def InitialCondition(ic, xyz):
for i in range(0,xyz.size1()):
x = xyz.get(i,0)
y = xyz.get(i,1)
z = xyz.get(i,2)
val=exp(-(x*x+y*y)*2) * exp(-((z+0.25)*(z+0.25)) * 50)
ic.set(i,0,val)
def InitialCondition(cm):
X = cm.get(xyz)
x = X[:,[0]]
y = X[:,[1]]
z = X[:,[2]]
return numpy.exp(-(x*x+y*y)*2) * numpy.exp(-((z+0.25)*(z+0.25)) * 50)
"""
Function for Advection parameters
"""
def Velocity(vel, xyz):
for i in range(0,xyz.size1()):
x = xyz.get(i,0)
y = xyz.get(i,1)
vel.set(i,0,Vx)
vel.set(i,1,Vy)
vel.set(i,2,Vz)
def Velocity(cm):
v = numpy.ndarray((cm.nEvaluationPoint(), 3))
v[:, 0] = Vx
v[:, 1] = Vy
v[:, 2] = Vz
return v
print('*** Loading the Mesh ***')
mesh = dgMesh(FileName, 0, [], False, 'periodicMesh.txt')
......@@ -54,7 +53,7 @@ groups=dgGroupCollection(mesh, order)
print ('*** Define Conservation Law ***')
xyz = groups.getFunctionCoordinates()
vel = functionPython(3, Velocity, [xyz])