/[escript]/branches/doubleplusgood/dudley/src/Assemble_setNormal.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/dudley/src/Assemble_setNormal.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26