/[escript]/branches/domexper/dudley/src/Assemble_setNormal.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Assemble_setNormal.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3221 - (show annotations)
Wed Sep 29 01:00:21 2010 UTC (8 years, 6 months ago) by jfenwick
File MIME type: text/plain
File size: 3636 byte(s)
Comment stripping

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* assemblage routines: calculates the normal vector at quadrature points on face elements */
18
19 /**************************************************************/
20
21 #include "Assemble.h"
22 #include "Util.h"
23 #ifdef _OPENMP
24 #include <omp.h>
25 #endif
26
27 #include "ShapeTable.h"
28
29 /**************************************************************/
30
31
32 void Dudley_Assemble_setNormal(Dudley_NodeFile* nodes, Dudley_ElementFile* elements, escriptDataC* normal) {
33 double *local_X=NULL, *dVdv=NULL,*normal_array;
34 index_t sign;
35 dim_t e,q, NN, NS, numDim, numQuad, numDim_local;
36 bool_t reduced_integration;
37 const double* dSdv=0;
38 if (nodes==NULL || elements==NULL) return;
39
40 switch (elements->numDim)
41 {
42 case 2: dSdv=&(DTDV_2D_alt[0][0]);break;
43 case 3: dSdv=&(DTDV_3D[0][0]);break;
44 default:
45 dSdv=&(DTDV_1D[0][0]);break;
46 }
47 Dudley_resetError();
48 NN=elements->numNodes;
49 numDim=nodes->numDim;
50 reduced_integration = Dudley_Assemble_reducedIntegrationOrder(normal);
51 numQuad=(!reduced_integration)?(elements->numDim+1):1;
52 numDim_local=elements->numDim;
53 NS=elements->numDim+1;
54
55 /* set some parameters */
56
57 sign=1;
58 /* check the dimensions of normal */
59 if (! (numDim==numDim_local || numDim-1==numDim_local)) {
60 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_setNormal: Cannot calculate normal vector");
61 } else if (! isDataPointShapeEqual(normal,1,&(numDim))) {
62 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_setNormal: illegal number of samples of normal Data object");
63 } else if (! numSamplesEqual(normal,numQuad,elements->numElements)) {
64 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_setNormal: illegal number of samples of normal Data object");
65 } else if (! isDataPointShapeEqual(normal,1,&(numDim))) {
66 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_setNormal: illegal point data shape of normal Data object");
67 } else if (!isExpanded(normal)) {
68 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_setNormal: expanded Data object is expected for normal.");
69 }
70
71 /* now we can start */
72 if (Dudley_noError()) {
73 requireWrite(normal);
74 #pragma omp parallel private(local_X,dVdv)
75 {
76 local_X=dVdv=NULL;
77 /* allocation of work arrays */
78 local_X=THREAD_MEMALLOC(NS*numDim,double);
79 dVdv=THREAD_MEMALLOC(numQuad*numDim*numDim_local,double);
80 if (!(Dudley_checkPtr(local_X) || Dudley_checkPtr(dVdv) ) ) {
81 /* open the element loop */
82 #pragma omp for private(e,q,normal_array) schedule(static)
83 for(e=0;e<elements->numElements;e++) {
84 /* gather local coordinates of nodes into local_X: */
85 Dudley_Util_Gather_double(NS,&(elements->Nodes[INDEX2(0,e,NN)]),numDim,nodes->Coordinates,local_X);
86 /* calculate dVdv(i,j,q)=local_X(i,n)*DSDv(n,j,q) */
87 Dudley_Util_SmallMatMult(numDim,numDim_local*numQuad,dVdv,NS,local_X,/*reference_element->BasisFunctions->*/dSdv);
88 /* get normalized vector: */
89 normal_array=getSampleDataRW(normal,e);
90 Dudley_NormalVector(numQuad,numDim,numDim_local,dVdv,normal_array);
91 for (q=0;q<numQuad*numDim;q++) normal_array[q]*=sign;
92 }
93 }
94 THREAD_MEMFREE(local_X);
95 THREAD_MEMFREE(dVdv);
96 }
97 }
98 }
99

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26