| 1 |
/* $Id$ */ |
| 2 |
|
| 3 |
/**************************************************************/ |
| 4 |
|
| 5 |
/* assemblage routines: calculate the gradient of nodal data at quadrature points */ |
| 6 |
|
| 7 |
/**************************************************************/ |
| 8 |
|
| 9 |
/* Copyrights by ACcESS Australia, 2003,2004 */ |
| 10 |
/* author: gross@access.edu.au */ |
| 11 |
/* Version: $Id$ */ |
| 12 |
|
| 13 |
/**************************************************************/ |
| 14 |
|
| 15 |
#include "escript/Data/DataC.h" |
| 16 |
#include "Common.h" |
| 17 |
#include "Finley.h" |
| 18 |
#include "Assemble.h" |
| 19 |
#include "NodeFile.h" |
| 20 |
#include "ElementFile.h" |
| 21 |
#include "Util.h" |
| 22 |
#ifdef _OPENMP |
| 23 |
#include <omp.h> |
| 24 |
#endif |
| 25 |
/*****************************************************************/ |
| 26 |
|
| 27 |
|
| 28 |
void Finley_Assemble_gradient(Finley_NodeFile* nodes, Finley_ElementFile* elements, |
| 29 |
escriptDataC* grad_data,escriptDataC* data) { |
| 30 |
|
| 31 |
double *local_X=NULL, *local_data=NULL, *dVdv=NULL, *dvdV=NULL, *Vol=NULL, *d_datadv=NULL, *gradS=NULL,*data_array; |
| 32 |
int numNodes,e,node_offset,*resort_nodes,i,q,type,NS_DOF,NN_DOF,dof_offset; |
| 33 |
#define NODES 0 |
| 34 |
#define DOF 1 |
| 35 |
#define REDUCED_DOF 2 |
| 36 |
if (nodes==NULL || elements==NULL) return; |
| 37 |
int NN=elements->ReferenceElement->Type->numNodes; |
| 38 |
int NS=elements->ReferenceElement->Type->numShapes; |
| 39 |
int id[NN]; |
| 40 |
int numDim=nodes->numDim; |
| 41 |
int data_type=getFunctionSpaceType(data); |
| 42 |
int numComps=getDataPointSize(data); |
| 43 |
int numQuad=elements->ReferenceElement->numQuadNodes; |
| 44 |
for (i=0;i<NN;i++) id[i]=i; |
| 45 |
|
| 46 |
/* set some parameter */ |
| 47 |
|
| 48 |
if (data_type==FINLEY_NODES) { |
| 49 |
type=NODES; |
| 50 |
resort_nodes=id; |
| 51 |
NN_DOF=elements->ReferenceElement->Type->numNodes; |
| 52 |
NS_DOF=elements->ReferenceElement->Type->numShapes; |
| 53 |
gradS=elements->ReferenceElement->dSdv; |
| 54 |
numNodes=nodes->numNodes; |
| 55 |
} else if (data_type==FINLEY_DEGREES_OF_FREEDOM) { |
| 56 |
type=DOF; |
| 57 |
resort_nodes=id; |
| 58 |
NN_DOF=elements->ReferenceElement->Type->numNodes; |
| 59 |
NS_DOF=elements->ReferenceElement->Type->numShapes; |
| 60 |
gradS=elements->ReferenceElement->dSdv; |
| 61 |
numNodes=nodes->numDegreesOfFreedom; |
| 62 |
} else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) { |
| 63 |
type=REDUCED_DOF; |
| 64 |
resort_nodes=elements->ReferenceElement->Type->linearNodes; |
| 65 |
NN_DOF=elements->LinearReferenceElement->Type->numNodes; |
| 66 |
NS_DOF=elements->LinearReferenceElement->Type->numShapes; |
| 67 |
gradS=elements->LinearReferenceElement->dSdv; |
| 68 |
numNodes=nodes->reducedNumDegreesOfFreedom; |
| 69 |
} else { |
| 70 |
Finley_ErrorCode=TYPE_ERROR; |
| 71 |
sprintf(Finley_ErrorMsg,"Cannot calculate gradient of data"); |
| 72 |
} |
| 73 |
if (getFunctionSpaceType(grad_data)==FINLEY_CONTACT_ELEMENTS_2) { |
| 74 |
node_offset=NN-NS; |
| 75 |
dof_offset=NN_DOF-NS_DOF; |
| 76 |
} else { |
| 77 |
node_offset=0; |
| 78 |
dof_offset=0; |
| 79 |
} |
| 80 |
|
| 81 |
/* check the dimensions of interpolated_data and data */ |
| 82 |
|
| 83 |
if (numDim!=elements->ReferenceElement->Type->numDim) { |
| 84 |
Finley_ErrorCode=TYPE_ERROR; |
| 85 |
sprintf(Finley_ErrorMsg,"Gradient: Spatial and element dimension must match."); |
| 86 |
} else if (! numSamplesEqual(grad_data,numQuad,elements->numElements)) { |
| 87 |
Finley_ErrorCode=TYPE_ERROR; |
| 88 |
sprintf(Finley_ErrorMsg,"illegal number of samples in gradient Data object"); |
| 89 |
} else if (! numSamplesEqual(data,1,numNodes)) { |
| 90 |
Finley_ErrorCode=TYPE_ERROR; |
| 91 |
sprintf(Finley_ErrorMsg,"illegal number of samples of input Data object"); |
| 92 |
} else if (numDim*numComps!=getDataPointSize(grad_data)) { |
| 93 |
Finley_ErrorCode=TYPE_ERROR; |
| 94 |
sprintf(Finley_ErrorMsg,"illegal number of components in gradient data object."); |
| 95 |
} else if (!isExpanded(grad_data)) { |
| 96 |
Finley_ErrorCode=TYPE_ERROR; |
| 97 |
sprintf(Finley_ErrorMsg,"expanded Data object is expected for output data."); |
| 98 |
} |
| 99 |
|
| 100 |
/* now we can start */ |
| 101 |
|
| 102 |
if (Finley_ErrorCode==NO_ERROR) { |
| 103 |
#pragma omp parallel private(local_X,local_data,dvdV,dVdv,Vol,d_datadv) |
| 104 |
{ |
| 105 |
local_X=local_data=dVdv=dvdV=Vol=d_datadv=NULL; |
| 106 |
/* allocation of work arrays */ |
| 107 |
local_X=THREAD_MEMALLOC(NS*numDim,double); |
| 108 |
local_data=THREAD_MEMALLOC(NS*numComps,double); |
| 109 |
dVdv=THREAD_MEMALLOC(numQuad*numDim*numDim,double); |
| 110 |
dvdV=THREAD_MEMALLOC(numQuad*numDim*numDim,double); |
| 111 |
Vol=THREAD_MEMALLOC(numQuad,double); |
| 112 |
d_datadv=THREAD_MEMALLOC(numQuad*numComps*numDim,double); |
| 113 |
if (!(Finley_checkPtr(local_X) || Finley_checkPtr(dVdv) || Finley_checkPtr(dvdV) || Finley_checkPtr(Vol) || Finley_checkPtr(d_datadv) || Finley_checkPtr(local_data) )) { |
| 114 |
/* open the element loop */ |
| 115 |
#pragma omp for private(e,i,q,data_array) schedule(static) |
| 116 |
for(e=0;e<elements->numElements;e++) { |
| 117 |
/* gather local coordinates of nodes into local_X: */ |
| 118 |
Finley_Util_Gather_double(NS,&(elements->Nodes[INDEX2(node_offset,e,NN)]),numDim,nodes->Coordinates,local_X); |
| 119 |
/* calculate dVdv(i,j,q)=local_X(i,n)*DSDv(n,j,q) */ |
| 120 |
Finley_Util_SmallMatMult(numDim,numDim*numQuad,dVdv,NS,local_X,elements->ReferenceElement->dSdv); |
| 121 |
/* dvdV=invert(dVdv) */ |
| 122 |
Finley_Util_InvertSmallMat(numQuad,numDim,dVdv,dvdV,Vol); |
| 123 |
/* gather local data into local_data(numComps,NS_DOF): */ |
| 124 |
switch (type) { |
| 125 |
case NODES: |
| 126 |
for (q=0;q<NS_DOF;q++) { |
| 127 |
i=elements->Nodes[INDEX2(resort_nodes[dof_offset+q],e,NN)]; |
| 128 |
data_array=getSampleData(data,i); |
| 129 |
Finley_copyDouble(numComps,data_array,local_data+q*numComps); |
| 130 |
} |
| 131 |
break; |
| 132 |
case DOF: |
| 133 |
for (q=0;q<NS_DOF;q++) { |
| 134 |
i=elements->Nodes[INDEX2(resort_nodes[dof_offset+q],e,NN)]; |
| 135 |
data_array=getSampleData(data,nodes->degreeOfFreedom[i]); |
| 136 |
Finley_copyDouble(numComps,data_array,local_data+q*numComps); |
| 137 |
|
| 138 |
} |
| 139 |
break; |
| 140 |
case REDUCED_DOF: |
| 141 |
for (q=0;q<NS_DOF;q++) { |
| 142 |
i=elements->Nodes[INDEX2(resort_nodes[dof_offset+q],e,NN)]; |
| 143 |
data_array=getSampleData(data,nodes->reducedDegreeOfFreedom[i]); |
| 144 |
Finley_copyDouble(numComps,data_array,local_data+q*numComps); |
| 145 |
} |
| 146 |
break; |
| 147 |
} |
| 148 |
/* calculate d_datadv(l,i,q)=local_data(l,n)*DSDv(n,i,q) */ |
| 149 |
Finley_Util_SmallMatMult(numComps,numDim*numQuad,d_datadv,NS_DOF,local_data,gradS); |
| 150 |
/* calculate grad_data(l,i)=d_datadv(l,i,q)*dvdV(i,k,q) */ |
| 151 |
Finley_Util_SmallMatSetMult(numQuad,numComps,numDim,getSampleData(grad_data,e),numDim,d_datadv,dvdV); |
| 152 |
} /* for */ |
| 153 |
} |
| 154 |
THREAD_MEMFREE(local_X); |
| 155 |
THREAD_MEMFREE(dVdv); |
| 156 |
THREAD_MEMFREE(dvdV); |
| 157 |
THREAD_MEMFREE(Vol); |
| 158 |
THREAD_MEMFREE(local_data); |
| 159 |
THREAD_MEMFREE(d_datadv); |
| 160 |
} |
| 161 |
} |
| 162 |
#undef NODES |
| 163 |
#undef DOF |
| 164 |
#undef REDUCED_DOF |
| 165 |
} |
| 166 |
/* |
| 167 |
* $Log$ |
| 168 |
* Revision 1.4 2004/12/15 07:08:32 jgs |
| 169 |
* *** empty log message *** |
| 170 |
* |
| 171 |
* |
| 172 |
* |
| 173 |
*/ |
Parent Directory
| 
Revision Log