finley/src/Assemble_getSize.c


/*******************************************************
*
* Copyright (c) 2003-2009 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/


/**************************************************************/

/*    assemblage routines: */

/*    calculates the minimum distance between two vertices of elements and assigns the value to each  */
/*    quadrature point in element_size                                                                         */


/**************************************************************/

#include "Assemble.h"
#include "Util.h"
#ifdef _OPENMP
#include <omp.h>
#endif

/**************************************************************/
void Finley_Assemble_getSize(Finley_NodeFile* nodes, Finley_ElementFile* elements, escriptDataC* element_size) {

  Finley_ShapeFunction *shape=NULL;
  Finley_ReferenceElement*  refElement=NULL;
  double *local_X=NULL,*element_size_array;
  dim_t e,n0,n1,q,i, NVertices, NN, NS, numQuad, numDim;
  index_t node_offset;
  double d,diff,min_diff, f;
  Finley_resetError();

  if (nodes==NULL || elements==NULL) return;
  refElement=Finley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet,Finley_Assemble_reducedIntegrationOrder(element_size));
  shape= refElement->Parametrization;
  numDim=nodes->numDim;
  numQuad=shape->numQuadNodes;
  NN=elements->numNodes;
  NS=refElement->Parametrization->Type->numShapes;
  NVertices=refElement->Parametrization->Type->numVertices;
  
  
  /* set a few more parameters */

  if (getFunctionSpaceType(element_size)==FINLEY_CONTACT_ELEMENTS_2) {
      node_offset=refElement->Type->offsets[1];
  } else {
      node_offset=refElement->Type->offsets[0];
  }
  f=pow(0.5, pow((double)(refElement->Type->numSubElements), 1./(double)(numDim))-1);

  /* check the dimensions of element_size */

  if (! numSamplesEqual(element_size,numQuad,elements->numElements)) {
       Finley_setError(TYPE_ERROR,"Finley_Assemble_getSize: illegal number of samples of element size Data object");
  } else if (! isDataPointShapeEqual(element_size,0,&(numDim))) {
       Finley_setError(TYPE_ERROR,"Finley_Assemble_getSize: illegal data point shape of element size Data object");
  }  else if (!isExpanded(element_size)) {
       Finley_setError(TYPE_ERROR,"Finley_Assemble_getSize: expanded Data object is expected for element size.");
  }
  /* now we can start: */

  if (Finley_noError()) {
    requireWrite(element_size);
        #pragma omp parallel private(local_X)
        {
            /* allocation of work arrays */
            local_X=THREAD_MEMALLOC(NN*numDim,double);
            if (! Finley_checkPtr(local_X) ) {
                /* open the element loop */
                #pragma omp for private(e,min_diff,diff,n0,n1,d,q,i,element_size_array) schedule(static)
                for(e=0;e<elements->numElements;e++) {
                    /* gather local coordinates of nodes into local_X(numDim,NN): */
                    Finley_Util_Gather_double(NS,&(elements->Nodes[INDEX2(node_offset,e,NN)]),numDim,nodes->Coordinates,local_X);
                    /* calculate minimal differences */
                    min_diff=-1;
                    for (n0=0;n0<NVertices;n0++) {
                        for (n1=n0+1;n1<NVertices;n1++) {
                            diff=0;
                            for (i=0;i<numDim;i++) {
                                d=local_X[INDEX2(i,n0,numDim)]-local_X[INDEX2(i,n1,numDim)];
                                diff+=d*d;
                            }
                            if (min_diff<0) {
                                min_diff=diff;
                            } else {
                                min_diff=MIN(min_diff,diff);
                            }
                        }
                    }
                    min_diff=sqrt(MAX(min_diff,0))*f;
                    /* set all values to min_diff */
                    element_size_array=getSampleDataRW(element_size,e);
                    for (q=0;q<numQuad;q++) element_size_array[q]=min_diff;
                }
            }
            THREAD_MEMFREE(local_X);
       } /* end of parallel region */
  }
  return;
}
1
2	/*******************************************************
3	*
4	* Copyright (c) 2003-2009 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: */
18
19	/* calculates the minimum distance between two vertices of elements and assigns the value to each */
20	/* quadrature point in element_size */
21
22
23	/**************************************************************/
24
25	#include "Assemble.h"
26	#include "Util.h"
27	#ifdef _OPENMP
28	#include <omp.h>
29	#endif
30
31	/**************************************************************/
32	void Finley_Assemble_getSize(Finley_NodeFile* nodes, Finley_ElementFile* elements, escriptDataC* element_size) {
33
34	Finley_ShapeFunction *shape=NULL;
35	Finley_ReferenceElement* refElement=NULL;
36	double local_X=NULL,element_size_array;
37	dim_t e,n0,n1,q,i, NVertices, NN, NS, numQuad, numDim;
38	index_t node_offset;
39	double d,diff,min_diff, f;
40	Finley_resetError();
41
42	if (nodes==NULL \|\| elements==NULL) return;
43	refElement=Finley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet,Finley_Assemble_reducedIntegrationOrder(element_size));
44	shape= refElement->Parametrization;
45	numDim=nodes->numDim;
46	numQuad=shape->numQuadNodes;
47	NN=elements->numNodes;
48	NS=refElement->Parametrization->Type->numShapes;
49	NVertices=refElement->Parametrization->Type->numVertices;
50
51
52	/* set a few more parameters */
53
54	if (getFunctionSpaceType(element_size)==FINLEY_CONTACT_ELEMENTS_2) {
55	node_offset=refElement->Type->offsets[1];
56	} else {
57	node_offset=refElement->Type->offsets[0];
58	}
59	f=pow(0.5, pow((double)(refElement->Type->numSubElements), 1./(double)(numDim))-1);
60
61	/* check the dimensions of element_size */
62
63	if (! numSamplesEqual(element_size,numQuad,elements->numElements)) {
64	Finley_setError(TYPE_ERROR,"Finley_Assemble_getSize: illegal number of samples of element size Data object");
65	} else if (! isDataPointShapeEqual(element_size,0,&(numDim))) {
66	Finley_setError(TYPE_ERROR,"Finley_Assemble_getSize: illegal data point shape of element size Data object");
67	} else if (!isExpanded(element_size)) {
68	Finley_setError(TYPE_ERROR,"Finley_Assemble_getSize: expanded Data object is expected for element size.");
69	}
70	/* now we can start: */
71
72	if (Finley_noError()) {
73	requireWrite(element_size);
74	#pragma omp parallel private(local_X)
75	{
76	/* allocation of work arrays */
77	local_X=THREAD_MEMALLOC(NN*numDim,double);
78	if (! Finley_checkPtr(local_X) ) {
79	/* open the element loop */
80	#pragma omp for private(e,min_diff,diff,n0,n1,d,q,i,element_size_array) schedule(static)
81	for(e=0;e<elements->numElements;e++) {
82	/* gather local coordinates of nodes into local_X(numDim,NN): */
83	Finley_Util_Gather_double(NS,&(elements->Nodes[INDEX2(node_offset,e,NN)]),numDim,nodes->Coordinates,local_X);
84	/* calculate minimal differences */
85	min_diff=-1;
86	for (n0=0;n0<NVertices;n0++) {
87	for (n1=n0+1;n1<NVertices;n1++) {
88	diff=0;
89	for (i=0;i<numDim;i++) {
90	d=local_X[INDEX2(i,n0,numDim)]-local_X[INDEX2(i,n1,numDim)];
91	diff+=d*d;
92	}
93	if (min_diff<0) {
94	min_diff=diff;
95	} else {
96	min_diff=MIN(min_diff,diff);
97	}
98	}
99	}
100	min_diff=sqrt(MAX(min_diff,0))*f;
101	/* set all values to min_diff */
102	element_size_array=getSampleDataRW(element_size,e);
103	for (q=0;q<numQuad;q++) element_size_array[q]=min_diff;
104	}
105	}
106	THREAD_MEMFREE(local_X);
107	} /* end of parallel region */
108	}
109	return;
110	}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision