From 758ed62c1512d9e731fa12d2774577bef885bf9e Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Fri, 18 Jun 2010 12:17:55 +0000
Subject: [PATCH] bindings in polynomialBasis and PViewData
---
Common/LuaBindings.cpp | 2 ++
Geo/MLine.cpp | 20 ++++++-------
Geo/MPoint.h | 2 +-
Geo/MPrism.cpp | 8 ++---
Geo/MQuadrangle.cpp | 40 ++++++++++++-------------
Geo/MTetrahedron.cpp | 30 +++++++++----------
Geo/MTriangle.cpp | 40 ++++++++++++-------------
Mesh/HighOrder.cpp | 10 +++----
Numeric/polynomialBasis.cpp | 22 ++++++++++----
Numeric/polynomialBasis.h | 15 +++++++++-
Post/PViewData.cpp | 53 ++++++++++++++++++++++++++++++++-
Post/PViewData.h | 4 +++
Post/PViewDataList.cpp | 2 +-
benchmarks/elasticity/conge.lua | 38 ++++++++++++++++++-----
14 files changed, 196 insertions(+), 90 deletions(-)
diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 89c4f644be..82910aaa5f 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -30,6 +30,7 @@
#include "GmshMessage.h"
#include "linearSystem.h"
#include "Options.h"
+#include "polynomialBasis.h"
#if defined(HAVE_OPENGL)
#include "drawContext.h"
@@ -412,6 +413,7 @@ binding::binding()
gmshOptions::registerBindings(this);
Msg::registerBindings(this);
linearSystem<double>::registerBindings(this);
+ polynomialBasis::registerBindings(this);
#if defined(HAVE_SOLVER)
function::registerBindings(this);
linearSystemCSRGmm<double>::registerBindings(this);
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index efdd93e4b3..36e1bb5cdf 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -14,16 +14,16 @@ const polynomialBasis* MLine::getFunctionSpace(int o) const
int order = (o == -1) ? getPolynomialOrder() : o;
switch (order) {
- case 1: return &polynomialBases::find(MSH_LIN_2);
- case 2: return &polynomialBases::find(MSH_LIN_3);
- case 3: return &polynomialBases::find(MSH_LIN_4);
- case 4: return &polynomialBases::find(MSH_LIN_5);
- case 5: return &polynomialBases::find(MSH_LIN_6);
- case 6: return &polynomialBases::find(MSH_LIN_7);
- case 7: return &polynomialBases::find(MSH_LIN_8);
- case 8: return &polynomialBases::find(MSH_LIN_9);
- case 9: return &polynomialBases::find(MSH_LIN_10);
- case 10: return &polynomialBases::find(MSH_LIN_11);
+ case 1: return polynomialBases::find(MSH_LIN_2);
+ case 2: return polynomialBases::find(MSH_LIN_3);
+ case 3: return polynomialBases::find(MSH_LIN_4);
+ case 4: return polynomialBases::find(MSH_LIN_5);
+ case 5: return polynomialBases::find(MSH_LIN_6);
+ case 6: return polynomialBases::find(MSH_LIN_7);
+ case 7: return polynomialBases::find(MSH_LIN_8);
+ case 8: return polynomialBases::find(MSH_LIN_9);
+ case 9: return polynomialBases::find(MSH_LIN_10);
+ case 10: return polynomialBases::find(MSH_LIN_11);
default: Msg::Error("Order %d line function space not implemented", order);
}
return 0;
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index 75c6183938..81b6f2cb3d 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -51,7 +51,7 @@ class MPoint : public MElement {
s[0][0] = s[0][1] = s[0][2] = 0.;
}
- virtual const polynomialBasis* getFunctionSpace(int o) const { return &polynomialBases::find(MSH_PNT); }
+ virtual const polynomialBasis* getFunctionSpace(int o) const { return polynomialBases::find(MSH_PNT); }
virtual bool isInside(double u, double v, double w)
{
return true;
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index e245eaeefc..d6e99b7063 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -38,15 +38,15 @@ const polynomialBasis* MPrism::getFunctionSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 1: return &polynomialBases::find(MSH_PRI_6);
- case 2: return &polynomialBases::find(MSH_PRI_18);
+ case 1: return polynomialBases::find(MSH_PRI_6);
+ case 2: return polynomialBases::find(MSH_PRI_18);
default: Msg::Error("Order %d prism function space not implemented", order);
}
}
else {
switch (order) {
- case 1: return &polynomialBases::find(MSH_PRI_6);
- case 2: return &polynomialBases::find(MSH_PRI_18);
+ case 1: return polynomialBases::find(MSH_PRI_6);
+ case 2: return polynomialBases::find(MSH_PRI_18);
default: Msg::Error("Order %d prism function space not implemented", order);
}
}
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index e5d4afe433..92eacb3a1b 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -24,29 +24,29 @@ const polynomialBasis* MQuadrangle::getFunctionSpace(int o) const
if ((nf == 0) && (o == -1)) {
switch (order) {
- case 1: return &polynomialBases::find(MSH_QUA_4);
- case 2: return &polynomialBases::find(MSH_QUA_8);
- case 3: return &polynomialBases::find(MSH_QUA_12);
- case 4: return &polynomialBases::find(MSH_QUA_16I);
- case 5: return &polynomialBases::find(MSH_QUA_20);
- case 6: return &polynomialBases::find(MSH_QUA_24);
- case 7: return &polynomialBases::find(MSH_QUA_28);
- case 8: return &polynomialBases::find(MSH_QUA_32);
- case 9: return &polynomialBases::find(MSH_QUA_36I);
- case 10: return &polynomialBases::find(MSH_QUA_40);
+ case 1: return polynomialBases::find(MSH_QUA_4);
+ case 2: return polynomialBases::find(MSH_QUA_8);
+ case 3: return polynomialBases::find(MSH_QUA_12);
+ case 4: return polynomialBases::find(MSH_QUA_16I);
+ case 5: return polynomialBases::find(MSH_QUA_20);
+ case 6: return polynomialBases::find(MSH_QUA_24);
+ case 7: return polynomialBases::find(MSH_QUA_28);
+ case 8: return polynomialBases::find(MSH_QUA_32);
+ case 9: return polynomialBases::find(MSH_QUA_36I);
+ case 10: return polynomialBases::find(MSH_QUA_40);
}
}
switch (order) {
- case 1: return &polynomialBases::find(MSH_QUA_4);
- case 2: return &polynomialBases::find(MSH_QUA_9);
- case 3: return &polynomialBases::find(MSH_QUA_16);
- case 4: return &polynomialBases::find(MSH_QUA_25);
- case 5: return &polynomialBases::find(MSH_QUA_36);
- case 6: return &polynomialBases::find(MSH_QUA_49);
- case 7: return &polynomialBases::find(MSH_QUA_64);
- case 8: return &polynomialBases::find(MSH_QUA_81);
- case 9: return &polynomialBases::find(MSH_QUA_100);
- case 10: return &polynomialBases::find(MSH_QUA_121);
+ case 1: return polynomialBases::find(MSH_QUA_4);
+ case 2: return polynomialBases::find(MSH_QUA_9);
+ case 3: return polynomialBases::find(MSH_QUA_16);
+ case 4: return polynomialBases::find(MSH_QUA_25);
+ case 5: return polynomialBases::find(MSH_QUA_36);
+ case 6: return polynomialBases::find(MSH_QUA_49);
+ case 7: return polynomialBases::find(MSH_QUA_64);
+ case 8: return polynomialBases::find(MSH_QUA_81);
+ case 9: return polynomialBases::find(MSH_QUA_100);
+ case 10: return polynomialBases::find(MSH_QUA_121);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
return 0;
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index ca1cc08c55..2cdeb6704e 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -111,26 +111,26 @@ const polynomialBasis* MTetrahedron::getFunctionSpace(int o) const
if ((nv == 0) && (o == -1)) {
switch (order) {
- case 1: return &polynomialBases::find(MSH_TET_4);
- case 2: return &polynomialBases::find(MSH_TET_10);
- case 3: return &polynomialBases::find(MSH_TET_20);
- case 4: return &polynomialBases::find(MSH_TET_34);
- case 5: return &polynomialBases::find(MSH_TET_52);
+ case 1: return polynomialBases::find(MSH_TET_4);
+ case 2: return polynomialBases::find(MSH_TET_10);
+ case 3: return polynomialBases::find(MSH_TET_20);
+ case 4: return polynomialBases::find(MSH_TET_34);
+ case 5: return polynomialBases::find(MSH_TET_52);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
else {
switch (order) {
- case 1: return &polynomialBases::find(MSH_TET_4);
- case 2: return &polynomialBases::find(MSH_TET_10);
- case 3: return &polynomialBases::find(MSH_TET_20);
- case 4: return &polynomialBases::find(MSH_TET_35);
- case 5: return &polynomialBases::find(MSH_TET_56);
- case 6: return &polynomialBases::find(MSH_TET_84);
- case 7: return &polynomialBases::find(MSH_TET_120);
- case 8: return &polynomialBases::find(MSH_TET_165);
- case 9: return &polynomialBases::find(MSH_TET_220);
- case 10: return &polynomialBases::find(MSH_TET_286);
+ case 1: return polynomialBases::find(MSH_TET_4);
+ case 2: return polynomialBases::find(MSH_TET_10);
+ case 3: return polynomialBases::find(MSH_TET_20);
+ case 4: return polynomialBases::find(MSH_TET_35);
+ case 5: return polynomialBases::find(MSH_TET_56);
+ case 6: return polynomialBases::find(MSH_TET_84);
+ case 7: return polynomialBases::find(MSH_TET_120);
+ case 8: return polynomialBases::find(MSH_TET_165);
+ case 9: return polynomialBases::find(MSH_TET_220);
+ case 10: return polynomialBases::find(MSH_TET_286);
default: Msg::Error("Order %d tetrahedron function space not implemented", order);
}
}
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index ce58c06cec..78fc812fa1 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -73,31 +73,31 @@ const polynomialBasis* MTriangle::getFunctionSpace(int o) const
if ((nf == 0) && (o == -1)) {
switch (order) {
- case 1: return &polynomialBases::find(MSH_TRI_3);
- case 2: return &polynomialBases::find(MSH_TRI_6);
- case 3: return &polynomialBases::find(MSH_TRI_9);
- case 4: return &polynomialBases::find(MSH_TRI_12);
- case 5: return &polynomialBases::find(MSH_TRI_15I);
- case 6: return &polynomialBases::find(MSH_TRI_18);
- case 7: return &polynomialBases::find(MSH_TRI_21I);
- case 8: return &polynomialBases::find(MSH_TRI_24);
- case 9: return &polynomialBases::find(MSH_TRI_27);
- case 10: return &polynomialBases::find(MSH_TRI_30);
+ case 1: return polynomialBases::find(MSH_TRI_3);
+ case 2: return polynomialBases::find(MSH_TRI_6);
+ case 3: return polynomialBases::find(MSH_TRI_9);
+ case 4: return polynomialBases::find(MSH_TRI_12);
+ case 5: return polynomialBases::find(MSH_TRI_15I);
+ case 6: return polynomialBases::find(MSH_TRI_18);
+ case 7: return polynomialBases::find(MSH_TRI_21I);
+ case 8: return polynomialBases::find(MSH_TRI_24);
+ case 9: return polynomialBases::find(MSH_TRI_27);
+ case 10: return polynomialBases::find(MSH_TRI_30);
default: Msg::Error("Order %d triangle incomplete function space not implemented", order);
}
}
else {
switch (order) {
- case 1: return &polynomialBases::find(MSH_TRI_3);
- case 2: return &polynomialBases::find(MSH_TRI_6);
- case 3: return &polynomialBases::find(MSH_TRI_10);
- case 4: return &polynomialBases::find(MSH_TRI_15);
- case 5: return &polynomialBases::find(MSH_TRI_21);
- case 6: return &polynomialBases::find(MSH_TRI_28);
- case 7: return &polynomialBases::find(MSH_TRI_36);
- case 8: return &polynomialBases::find(MSH_TRI_45);
- case 9: return &polynomialBases::find(MSH_TRI_55);
- case 10: return &polynomialBases::find(MSH_TRI_66);
+ case 1: return polynomialBases::find(MSH_TRI_3);
+ case 2: return polynomialBases::find(MSH_TRI_6);
+ case 3: return polynomialBases::find(MSH_TRI_10);
+ case 4: return polynomialBases::find(MSH_TRI_15);
+ case 5: return polynomialBases::find(MSH_TRI_21);
+ case 6: return polynomialBases::find(MSH_TRI_28);
+ case 7: return polynomialBases::find(MSH_TRI_36);
+ case 8: return polynomialBases::find(MSH_TRI_45);
+ case 9: return polynomialBases::find(MSH_TRI_55);
+ case 10: return polynomialBases::find(MSH_TRI_66);
default: Msg::Error("Order %d triangle function space not implemented", order);
}
}
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index bdae7c5046..b6c5c82852 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -515,15 +515,15 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
switch (nPts){
case 2 :
- points = polynomialBases::find(MSH_TRI_10).points;
+ points = polynomialBases::find(MSH_TRI_10)->points;
start = 9;
break;
case 3 :
- points = polynomialBases::find(MSH_TRI_15).points;
+ points = polynomialBases::find(MSH_TRI_15)->points;
start = 12;
break;
case 4 :
- points = polynomialBases::find(MSH_TRI_21).points;
+ points = polynomialBases::find(MSH_TRI_21)->points;
start = 15;
break;
default :
@@ -608,11 +608,11 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
switch (nPts){
case 3 :
- points = polynomialBases::find(MSH_TET_35).points;
+ points = polynomialBases::find(MSH_TET_35)->points;
start = 34;
break;
case 4 :
- points = polynomialBases::find(MSH_TET_56).points;
+ points = polynomialBases::find(MSH_TET_56)->points;
start = 52;
break;
default :
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 79a94761f0..474f4f0250 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -906,10 +906,10 @@ static void generate1dVertexClosure(polynomialBasis::clCont &closure)
}
std::map<int, polynomialBasis> polynomialBases::fs;
-const polynomialBasis &polynomialBases::find(int tag)
+const polynomialBasis *polynomialBases::find(int tag)
{
std::map<int, polynomialBasis>::const_iterator it = fs.find(tag);
- if (it != fs.end()) return it->second;
+ if (it != fs.end()) return &it->second;
polynomialBasis F;
F.numFaces = -1;
@@ -1305,7 +1305,7 @@ const polynomialBasis &polynomialBases::find(int tag)
// }
fs.insert(std::make_pair(tag, F));
- return fs[tag];
+ return &fs[tag];
}
@@ -1317,8 +1317,8 @@ const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
std::map<std::pair<int, int>, fullMatrix<double> >::const_iterator it = injector.find(key);
if (it != injector.end()) return it->second;
- const polynomialBasis& fs1 = find(tag1);
- const polynomialBasis& fs2 = find(tag2);
+ const polynomialBasis& fs1 = *find(tag1);
+ const polynomialBasis& fs2 = *find(tag2);
fullMatrix<double> inj(fs1.points.size1(), fs2.points.size1());
@@ -1332,3 +1332,15 @@ const fullMatrix<double> &polynomialBases::findInjector(int tag1, int tag2)
injector.insert(std::make_pair(key, inj));
return injector[key];
}
+
+#include "Bindings.h"
+void polynomialBasis::registerBindings(binding *b) {
+ classBinding *cb = b->addClass<polynomialBasis>("polynomialBasis");
+ cb->setDescription("polynomial shape functions for elements");
+ methodBinding *mb = cb->addMethod("f",(void (polynomialBasis::*)(fullMatrix<double>&, fullMatrix<double>&))&polynomialBasis::f);
+ mb->setDescription("evaluate the shape functions");
+ mb->setArgNames("nodes","values",NULL);
+ mb = cb->addMethod("find",&polynomialBases::find);
+ mb->setDescription("return the polynomial basis corresponding to an element type");
+ mb->setArgNames("elementType",NULL);
+}
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 80cece6706..e0f38baa21 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -62,6 +62,7 @@ inline double pow_int (const double &a , const int &n) {
}
}
+class binding;
// presently those function spaces are only for simplices and quads;
// should be extended to other elements like hexes
@@ -102,6 +103,17 @@ class polynomialBasis
}
}
}
+ // I would favour this interface that allows optimizations (not implemented) and is easier to bind
+ inline void f(fullMatrix<double> &coord, fullMatrix<double> &sf) {
+ double p[256];
+ sf.resize (coefficients.size1(), coord.size1());
+ for (int iPoint=0; iPoint< coord.size1(); iPoint++) {
+ evaluateMonomials(coord(iPoint,0), coord(iPoint,1), coord(iPoint,2), p);
+ for (int i = 0; i < coefficients.size1(); i++)
+ for (int j = 0; j < coefficients.size2(); j++)
+ sf(i,iPoint) += coefficients(i, j) * p[j];
+ }
+ }
inline void df(double u, double v, double w, double grads[][3]) const
{
switch (monomials.size2()) {
@@ -242,6 +254,7 @@ class polynomialBasis
break;
}
}
+ static void registerBindings(binding *b);
};
class polynomialBases
@@ -250,7 +263,7 @@ class polynomialBases
static std::map<int, polynomialBasis> fs;
static std::map<std::pair<int, int>, fullMatrix<double> > injector;
public :
- static const polynomialBasis &find(int);
+ static const polynomialBasis *find(int);
static const fullMatrix<double> &findInjector(int, int);
};
diff --git a/Post/PViewData.cpp b/Post/PViewData.cpp
index 21da15f879..2a1dc767db 100644
--- a/Post/PViewData.cpp
+++ b/Post/PViewData.cpp
@@ -120,6 +120,27 @@ bool PViewData::combineSpace(nameData &nd)
return false;
}
+double PViewData::getValueBinding(int step, int ent, int ele, int nod, int comp)
+{
+ double val;
+ getValue(step,ent,ele,nod,comp,val);
+ return val;
+}
+
+void PViewData::getAllValuesForElementBinding(int step, int ent, int ele, fullMatrix<double> &m) {
+ int nNodes = getNumNodes(step,ent,ele);
+ int nComponents = getNumComponents(step,ent,ele);
+ for (int i=0; i<nNodes; i++)
+ for (int j=0; j<nComponents; j++)
+ getValue(step,ent,ele,i,j,m(i,j));
+}
+
+void PViewData::getAllNodesForElementBinding(int step, int ent, int ele, fullMatrix<double> &m) {
+ int nNodes = getNumNodes(step,ent,ele);
+ for (int i=0; i<nNodes; i++)
+ getNode(step,ent,ele,i,m(i,0),m(i,1),m(i,2));
+}
+
#include "Bindings.h"
void PViewData::registerBindings(binding *b) {
classBinding *cb = b->addClass<PViewData>("PViewData");
@@ -127,5 +148,35 @@ void PViewData::registerBindings(binding *b) {
methodBinding *cm;
cm = cb->addMethod("getNumEntities",&PViewData::getNumEntities);
cm->setArgNames("step",NULL);
- cm->setDescription("return the number of entities for a given time step");
+ cm->setDescription("return the number of entities for a given time step (-1 for default)");
+ cm = cb->addMethod("getNumElements",&PViewData::getNumElements);
+ cm->setArgNames("step","entity",NULL);
+ cm->setDescription("return the number of entities for a given time step (-1 for default) for a given entity (-1 for all)");
+ cm = cb->addMethod("getNumTriangles",&PViewData::getNumTriangles);
+ cm->setArgNames("step",NULL);
+ cm->setDescription("return the number of triangles for a given time step (-1 for default)");
+
+ cm = cb->addMethod("getNumNodes",&PViewData::getNumNodes);
+ cm->setArgNames("step","entity","element",NULL);
+ cm->setDescription("return the number of nodes of one element of an entity of a time step (-1 for default time step)");
+
+ cm = cb->addMethod("getNumValues",&PViewData::getNumValues);
+ cm->setArgNames("step","entity","element",NULL);
+ cm->setDescription("return the number of values of one element of an entity of a time step (-1 for default time step)");
+
+ cm = cb->addMethod("getNumComponents",&PViewData::getNumComponents);
+ cm->setArgNames("step","entity","element",NULL);
+ cm->setDescription("return the number of components of one element of an entity of a time step (-1 for default time step)");
+
+ cm = cb->addMethod("getValue",&PViewData::getValueBinding);
+ cm->setArgNames("step","entity","element","node","component",NULL);
+ cm->setDescription("return one of the values");
+
+ cm = cb->addMethod("getAllValuesForElement",&PViewData::getAllValuesForElementBinding);
+ cm->setArgNames("step","entity","element","values",NULL);
+ cm->setDescription("fill a fullMatrix with all values from the elements. The fullMatrix should have the size nbNodes x nbComponents.");
+
+ cm = cb->addMethod("getAllNodesForElement",&PViewData::getAllNodesForElementBinding);
+ cm->setArgNames("step","entity","element","coordinates",NULL);
+ cm->setDescription("fill a fullMatrix with all coordinates from the nodes of the elements. The fullMatrix should have the size nbNodes x 3");
}
diff --git a/Post/PViewData.h b/Post/PViewData.h
index 24879b7e10..5b0e24808f 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -132,6 +132,10 @@ class PViewData {
virtual void getValue(int step, int ent, int ele, int nod, int comp, double &val){}
virtual void setValue(int step, int ent, int ele, int nod, int comp, double val);
+ double getValueBinding(int step, int ent, int ele, int nod, int comp);
+ void getAllValuesForElementBinding(int step, int ent, int ele, fullMatrix<double> &m);
+ void getAllNodesForElementBinding(int step, int ent, int ele, fullMatrix<double> &m);
+
// return a scalar value (same as value for scalars, norm for
// vectors, etc.) associated with the node-th node from the ele-th
// element in the ent-th entity
diff --git a/Post/PViewDataList.cpp b/Post/PViewDataList.cpp
index 8c39b8d7f9..53a68f5ca2 100644
--- a/Post/PViewDataList.cpp
+++ b/Post/PViewDataList.cpp
@@ -784,7 +784,7 @@ void PViewDataList::setOrder2(int type)
case TYPE_PRI: typeMSH = MSH_PRI_18; break;
case TYPE_PYR: typeMSH = MSH_PYR_14; break;
}
- const polynomialBasis *fs = &polynomialBases::find(typeMSH);
+ const polynomialBasis *fs = polynomialBases::find(typeMSH);
if(!fs){
Msg::Error("Could not find function space for element type %d", typeMSH);
return;
diff --git a/benchmarks/elasticity/conge.lua b/benchmarks/elasticity/conge.lua
index 041ae640b0..1319eb85f5 100644
--- a/benchmarks/elasticity/conge.lua
+++ b/benchmarks/elasticity/conge.lua
@@ -15,13 +15,37 @@ sys = linearSystemCSRTaucs()
e:assemble(sys)
sys:systemSolve()
view = e:buildVonMisesView("vonMises")
-print(view)
-view:write("vonMises.msh",5,false)
+
view = e:buildDisplacementView("displacement")
view:write("displacement.msh",5,false)
-view = e:buildLagrangeMultiplierView("lagrangeMultiplier")
-view:write("lagrangeMultiplier.msh",5,false)
-view = e:buildElasticEnergyView("elasticEnergy")
-view:write("elasticEnergy.msh",5,false)
+
viewData = view:getData()
-print(viewData:getNumEntities(-1))
+nEntities = viewData:getNumEntities(-1)
+v = fullMatrix(3,3)
+nodes = fullMatrix(3,3)
+for i=0,nEntities-1 do
+ nElements = viewData:getNumElements(-1,i)
+ for j=0,nElements-1 do
+ if(viewData:getNumNodes(0,i,j)==3) then
+ viewData:getAllValuesForElement(0,i,j,v)
+ viewData:getAllNodesForElement(0,i,j,nodes)
+ end
+ end
+end
+
+basis = polynomialBasis.find(2)
+points = fullMatrix(3,3)
+points:set(0,0,0);
+points:set(0,1,0);
+points:set(0,2,0);
+
+points:set(1,0,1);
+points:set(1,1,0);
+points:set(1,2,0);
+
+points:set(2,0,0);
+points:set(2,1,1);
+points:set(2,2,0);
+
+basis:f(points,v)
+v:print("f")
--
GitLab