/[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 3090 - (show annotations)
Wed Aug 11 00:51:28 2010 UTC (9 years, 6 months ago) by jfenwick
File MIME type: text/plain
File size: 3668 byte(s)
Removed:
DUDLEY_CONTACT_ELEMENTS_1
DUDLEY_CONTACT_ELEMENTS_2
DUDLEY_REDUCED_CONTACT_ELEMENTS_1
DUDLEY_REDUCED_CONTACT_ELEMENTS_2

escript tests now query if the domain supports contact elements
before trying to use them.


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 /**************************************************************/
28
29
30 void Dudley_Assemble_setNormal(Dudley_NodeFile* nodes, Dudley_ElementFile* elements, escriptDataC* normal) {
31 double *local_X=NULL, *dVdv=NULL,*normal_array;
32 index_t sign,node_offset;
33 Dudley_ReferenceElement* reference_element=NULL;
34 dim_t e,q, NN, NS, numDim, numQuad, numDim_local;
35 if (nodes==NULL || elements==NULL) return;
36 Dudley_resetError();
37 reference_element=Dudley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, Dudley_Assemble_reducedIntegrationOrder(normal));
38 NN=elements->numNodes;
39 numDim=nodes->numDim;
40
41 numQuad=reference_element->Parametrization->numQuadNodes;
42 numDim_local=reference_element->Parametrization->Type->numDim;
43 NS=reference_element->Parametrization->Type->numShapes;
44
45
46
47 /* set some parameter */
48
49 node_offset=reference_element->Type->offsets[0];
50 sign=1;
51 /* check the dimensions of normal */
52 if (! (numDim==numDim_local || numDim-1==numDim_local)) {
53 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_setNormal: Cannot calculate normal vector");
54 } else if (! isDataPointShapeEqual(normal,1,&(numDim))) {
55 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_setNormal: illegal number of samples of normal Data object");
56 } else if (! numSamplesEqual(normal,numQuad,elements->numElements)) {
57 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_setNormal: illegal number of samples of normal Data object");
58 } else if (! isDataPointShapeEqual(normal,1,&(numDim))) {
59 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_setNormal: illegal point data shape of normal Data object");
60 } else if (!isExpanded(normal)) {
61 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_setNormal: expanded Data object is expected for normal.");
62 }
63
64 /* now we can start */
65 if (Dudley_noError()) {
66 requireWrite(normal);
67 #pragma omp parallel private(local_X,dVdv)
68 {
69 local_X=dVdv=NULL;
70 /* allocation of work arrays */
71 local_X=THREAD_MEMALLOC(NS*numDim,double);
72 dVdv=THREAD_MEMALLOC(numQuad*numDim*numDim_local,double);
73 if (!(Dudley_checkPtr(local_X) || Dudley_checkPtr(dVdv) ) ) {
74 /* open the element loop */
75 #pragma omp for private(e,q,normal_array) schedule(static)
76 for(e=0;e<elements->numElements;e++) {
77 /* gather local coordinates of nodes into local_X: */
78 Dudley_Util_Gather_double(NS,&(elements->Nodes[INDEX2(node_offset,e,NN)]),numDim,nodes->Coordinates,local_X);
79 /* calculate dVdv(i,j,q)=local_X(i,n)*DSDv(n,j,q) */
80 Dudley_Util_SmallMatMult(numDim,numDim_local*numQuad,dVdv,NS,local_X,reference_element->Parametrization->dSdv);
81 /* get normalized vector: */
82 normal_array=getSampleDataRW(normal,e);
83 Dudley_NormalVector(numQuad,numDim,numDim_local,dVdv,normal_array);
84 for (q=0;q<numQuad*numDim;q++) normal_array[q]*=sign;
85 }
86 }
87 THREAD_MEMFREE(local_X);
88 THREAD_MEMFREE(dVdv);
89 }
90 }
91 }
92

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26