1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#include "mex.h"
#include "mat.h"
#include "math.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    
    double *xg,*yg,*zg,*FF;
    long dimP, i, j,elt;
    int count=0;
    mwSize mwDimG[2];
    
    if(nrhs!=4 || nlhs!=1) {
        mexErrMsgTxt("[FF]=GetFiniteElementShapeFunctions3D(elt,xg,yg,zg)");
    }
    
    elt=mxGetScalar(prhs[0]);
    xg=mxGetData(prhs[1]);
    yg=mxGetData(prhs[2]);
    zg=mxGetData(prhs[3]);
    
    
    dimP=mxGetDimensions(prhs[1])[0];
    
    mwDimG[0] = dimP;
    mwDimG[1] = 8;
    plhs[0] = mxCreateNumericArray(2, mwDimG, mxDOUBLE_CLASS, mxREAL);
    FF= mxGetData(plhs[0]);

    
    
    
    switch (elt){
        case 8 :
            
                for(i=0;i<dimP;i++) {
                    *(FF + count)=0.125*(1-(*(xg+i)))*(1-(*(yg+i)))*(1-(*(zg+i)));
                    count++;
                }
                for(i=0;i<dimP;i++) {
                    *(FF + count)=0.125*(1+(*(xg+i)))*(1-(*(yg+i)))*(1-(*(zg+i)));
                    count++;
                }
                for(i=0;i<dimP;i++) {
                    *(FF + count)=0.125*(1+(*(xg+i)))*(1+(*(yg+i)))*(1-(*(zg+i)));
                    count++;
                }
                for(i=0;i<dimP;i++) {
                    *(FF + count)=0.125*(1-(*(xg+i)))*(1+(*(yg+i)))*(1-(*(zg+i)));
                    count++;
                }
                for(i=0;i<dimP;i++) {
                    *(FF + count)=0.125*(1-(*(xg+i)))*(1-(*(yg+i)))*(1+(*(zg+i)));
                    count++;
                }
                for(i=0;i<dimP;i++) {
                    *(FF + count)=0.125*(1+(*(xg+i)))*(1-(*(yg+i)))*(1+(*(zg+i)));
                    count++;
                }
                for(i=0;i<dimP;i++) {
                    *(FF + count)=0.125*(1+(*(xg+i)))*(1+(*(yg+i)))*(1+(*(zg+i)));
                    count++;
                }
                for(i=0;i<dimP;i++) {
                    *(FF + count)=0.125*(1-(*(xg+i)))*(1+(*(yg+i)))*(1+(*(zg+i)));
                    count++;
                }
            
            break;
            
        case 4:
                for(i=0;i<dimP;i++) {
                    *(FF + count)=1-(*(xg+i))-(*(yg+i))-(*(zg+i));
                    count++;
                }
                for(i=0;i<dimP;i++) {
                    *(FF + count)=(*(xg+i));
                    count++;
                }
                for(i=0;i<dimP;i++) {
                    *(FF + count)=(*(yg+i));
                    count++;
                }
                for(i=0;i<dimP;i++) {
                    *(FF + count)=(*(zg+i));
                    count++;
                }
            break;
        default :
            break;
            
    }
    
    
}