/[escript]/trunk/dudley/src/Assemble_setNormal.c
ViewVC logotype

Contents of /trunk/dudley/src/Assemble_setNormal.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years, 1 month ago) by jfenwick
File MIME type: text/plain
File size: 3768 byte(s)
Merging dudley and scons updates from branches

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 /* assemblage routines: calculates the normal vector at quadrature points on face elements */
17
18 /**************************************************************/
19
20 #include "Assemble.h"
21 #include "Util.h"
22 #ifdef _OPENMP
23 #include <omp.h>
24 #endif
25
26 #include "ShapeTable.h"
27
28 /**************************************************************/
29
30 void Dudley_Assemble_setNormal(Dudley_NodeFile * nodes, Dudley_ElementFile * elements, escriptDataC * normal)
31 {
32 double *local_X = NULL, *dVdv = NULL, *normal_array;
33 index_t sign;
34 dim_t e, q, NN, NS, numDim, numQuad, numDim_local;
35 bool_t reduced_integration;
36 const double *dSdv = 0;
37 if (nodes == NULL || elements == NULL)
38 return;
39
40 switch (elements->numDim)
41 {
42 case 2:
43 dSdv = &(DTDV_2D[0][0]);
44 break;
45 case 3:
46 dSdv = &(DTDV_3D[0][0]);
47 break;
48 default:
49 dSdv = &(DTDV_1D[0][0]);
50 break;
51 }
52 Dudley_resetError();
53 NN = elements->numNodes;
54 numDim = nodes->numDim;
55 reduced_integration = Dudley_Assemble_reducedIntegrationOrder(normal);
56 numQuad = (!reduced_integration) ? (elements->numDim + 1) : 1;
57 numDim_local = elements->numLocalDim;
58 NS = elements->numDim + 1;
59
60 /* set some parameters */
61
62 sign = 1;
63 /* check the dimensions of normal */
64 if (!(numDim == numDim_local || numDim - 1 == numDim_local))
65 {
66 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_setNormal: Cannot calculate normal vector");
67 }
68 else if (!isDataPointShapeEqual(normal, 1, &(numDim)))
69 {
70 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_setNormal: illegal number of samples of normal Data object");
71 }
72 else if (!numSamplesEqual(normal, numQuad, elements->numElements))
73 {
74 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_setNormal: illegal number of samples of normal Data object");
75 }
76 else if (!isDataPointShapeEqual(normal, 1, &(numDim)))
77 {
78 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_setNormal: illegal point data shape of normal Data object");
79 }
80 else if (!isExpanded(normal))
81 {
82 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_setNormal: expanded Data object is expected for normal.");
83 }
84
85 /* now we can start */
86 if (Dudley_noError())
87 {
88 requireWrite(normal);
89 #pragma omp parallel private(local_X,dVdv)
90 {
91 local_X = dVdv = NULL;
92 /* allocation of work arrays */
93 local_X = THREAD_MEMALLOC(NS * numDim, double);
94 dVdv = THREAD_MEMALLOC(numQuad * numDim * numDim_local, double);
95 if (!(Dudley_checkPtr(local_X) || Dudley_checkPtr(dVdv)))
96 {
97 /* open the element loop */
98 #pragma omp for private(e,q,normal_array) schedule(static)
99 for (e = 0; e < elements->numElements; e++)
100 {
101 /* gather local coordinates of nodes into local_X: */
102 Dudley_Util_Gather_double(NS, &(elements->Nodes[INDEX2(0, e, NN)]), numDim, nodes->Coordinates,
103 local_X);
104 /* calculate dVdv(i,j,q)=local_X(i,n)*DSDv(n,j,q) */
105 Dudley_Util_SmallMatMult(numDim, numDim_local * numQuad, dVdv, NS, local_X, dSdv);
106 /* get normalized vector: */
107 normal_array = getSampleDataRW(normal, e);
108 Dudley_NormalVector(numQuad, numDim, numDim_local, dVdv, normal_array);
109 for (q = 0; q < numQuad * numDim; q++)
110 normal_array[q] *= sign;
111 }
112 }
113 THREAD_MEMFREE(local_X);
114 THREAD_MEMFREE(dVdv);
115 }
116 }
117 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26