Home Research Papers Posters Manuals Talks Download People Help Wiki Internal

libs/ck-libs/datatransfer/test_GenericElement.C

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include "GenericElement.h"

class fixedConcreteElement : public ConcreteElementNodeData {
public:
        enum{nPts=4};
        CPoint pts[nPts], data[nPts];
        fixedConcreteElement(
                CPoint a,CPoint b,CPoint c,CPoint d)
        {
                pts[0]=a; pts[1]=b; pts[2]=c; pts[3]=d;
                
                for (int i=0;i<nPts;i++) data[i]=pts[i];
        }
        
        void shift(const CVector &v) {
                for (int i=0;i<nPts;i++) pts[i]+=v;
        }
        
        CPoint getNodeLocation(int i) const {return pts[i];}
        const double *getNodeData(int i) const {return data[i];}
        
        CPoint getCenter(void) const {
                CPoint r(0);
                for (int i=0;i<nPts;i++) r+=pts[i];
                return r*(1.0/nPts);
        }
};

void die(const char *fmt,const char *file,int line,const char *why) {
        fprintf(stderr,fmt,file,line,why);
        abort();
}

#define check(expr) { bool b=(expr); \
        if (!b) die("Check failed: %s:%d> (%s) is false\n",__FILE__,__LINE__,#expr); \
  }

bool equalVectors(const CVector &l,const CVector &r) {
        for (int i=0;i<3;i++)
                if (fabs(l[i]-r[i])>1.0e-10)
                        return false;
        return true;
}

int main() {
        CVector geoShift(1.2,2.3,3.4);
        
        fixedConcreteElement r(CPoint(0,0,0),CPoint(1,0,0),CPoint(0,1,0),CPoint(0,0,1));
        fixedConcreteElement l(CPoint(0,0,0),CPoint(-1,0,0),CPoint(0,1,0),CPoint(0,0,1));
        l.shift(geoShift);
        r.shift(geoShift);
        
        GenericElement tet(4);
        CVector natc;
        
        // Check element centers:
        check(tet.element_contains_point(r.getCenter(),r,natc));
        check(tet.element_contains_point(l.getCenter(),l,natc));
        check(!tet.element_contains_point(r.getCenter(),l,natc));
        check(!tet.element_contains_point(l.getCenter(),r,natc));
        
        // Check interpolation on centers:
        CVector interp;
        check(tet.interpolate(3,r,r.getCenter(),interp));
        check(equalVectors(interp,r.getCenter()-geoShift));
        check(tet.interpolate(3,l,l.getCenter(),interp));
        check(equalVectors(interp,l.getCenter()-geoShift));
        
        // Check boundaries:
        for (int i=0;i<4;i++) {
                check(tet.interpolate(3,r,r.getNodeLocation(i),interp));
                check(equalVectors(interp,CVector(r.getNodeData(i))));
                
                double frac=0.1+0.2*i;
                CPoint middle(0,frac,0); middle+=geoShift;
                CVector interpL,interpR;
                check(tet.interpolate(3,l,middle,interpL));
                check(tet.interpolate(3,r,middle,interpR));
                check(equalVectors(interpL,interpR));
                CVector theoryValue=(1-frac)*CVector(r.getNodeData(0))+
                                       frac*CVector(r.getNodeData(2));
                check(equalVectors(interpL,theoryValue));
        }
        
        printf("All tests passed.\n");
        return 0;
}

---