finley/src/ElementFile.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
*
*******************************************************/


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

/*   Finley: ElementFile */

/*   allocates an element file to hold elements of type id and with integration order order. */
/*   use Finley_Mesh_allocElementTable to allocate the element table (Id,Nodes,Tag,Owner). */

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

#include "ElementFile.h"

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

Finley_ElementFile* Finley_ElementFile_alloc(ElementTypeId id, index_t order, index_t reduced_order, Paso_MPIInfo *MPIInfo)
{
  extern Finley_RefElementInfo Finley_RefElement_InfoList[];
  dim_t NQ, reduced_NQ;
  Finley_ElementFile *out;
  
  /*   get the number of quadrature nodes needed to achieve integration order order: */
  
  if (order<0) order=MAX(2*Finley_RefElement_InfoList[id].numOrder,0);
  if (reduced_order<0) reduced_order=MAX(2*(Finley_RefElement_InfoList[id].numOrder-1),0);
  NQ= Finley_RefElement_InfoList[id].getNumQuadNodes(order);
  reduced_NQ= Finley_RefElement_InfoList[id].getNumQuadNodes(reduced_order);
  if (! Finley_noError()) return NULL;
  
  /*  allocate the return value */
  
  out=MEMALLOC(1,Finley_ElementFile);
  if (Finley_checkPtr(out)) return NULL;
  out->order = order;
  out->reduced_order = reduced_order;
  out->ReferenceElement=NULL;
  out->LinearReferenceElement=NULL;
  out->ReferenceElementReducedOrder=NULL;
  out->LinearReferenceElementReducedOrder=NULL;
  out->numElements=0;
  out->Id=NULL;
  out->Nodes=NULL;
  out->Tag=NULL;
  out->Color=NULL;
  out->minColor=0;
  out->maxColor=-1;
  out->jacobeans=NULL;
  out->jacobeans_reducedQ=NULL;
  out->jacobeans_reducedS=NULL;
  out->jacobeans_reducedS_reducedQ=NULL;

  out->Owner=NULL;                
  out->numTagsInUse=0;
  out->tagsInUse=NULL;

  out->MPIInfo = Paso_MPIInfo_getReference( MPIInfo );

  /*  allocate the reference element: */
  
  out->ReferenceElement=Finley_RefElement_alloc(id,NQ);
  out->jacobeans=Finley_ElementFile_Jacobeans_alloc(out->ReferenceElement);
  out->ReferenceElementReducedOrder=Finley_RefElement_alloc(id,reduced_NQ);
  out->jacobeans_reducedQ=Finley_ElementFile_Jacobeans_alloc(out->ReferenceElementReducedOrder);

  out->LinearReferenceElement=Finley_RefElement_alloc(Finley_RefElement_InfoList[id].LinearTypeId,NQ);
  out->jacobeans_reducedS=Finley_ElementFile_Jacobeans_alloc(out->LinearReferenceElement);
  out->LinearReferenceElementReducedOrder=Finley_RefElement_alloc(Finley_RefElement_InfoList[id].LinearTypeId,reduced_NQ);
  out->jacobeans_reducedS_reducedQ=Finley_ElementFile_Jacobeans_alloc(out->LinearReferenceElementReducedOrder);

  out->numNodes=out->ReferenceElement->Type->numNodes;

  if (! Finley_noError()) {
     Finley_ElementFile_free(out);
     return NULL;
  }
  return out;
}

/*  deallocates an element file: */

void Finley_ElementFile_free(Finley_ElementFile* in) {
  if (in!=NULL) {
     #ifdef Finley_TRACE
     if (in->ReferenceElement!=NULL) printf("element file for %s is deallocated.\n",in->ReferenceElement->Type->Name);
     #endif
     Finley_ElementFile_freeTable(in);   
     Finley_RefElement_dealloc(in->ReferenceElement);
     Finley_RefElement_dealloc(in->ReferenceElementReducedOrder);
     Finley_RefElement_dealloc(in->LinearReferenceElement);
     Finley_RefElement_dealloc(in->LinearReferenceElementReducedOrder);
     Finley_ElementFile_Jacobeans_dealloc(in->jacobeans);
     Finley_ElementFile_Jacobeans_dealloc(in->jacobeans_reducedS);
     Finley_ElementFile_Jacobeans_dealloc(in->jacobeans_reducedQ);
     Finley_ElementFile_Jacobeans_dealloc(in->jacobeans_reducedS_reducedQ);
     Paso_MPIInfo_free( in->MPIInfo );
     MEMFREE(in);      
  }
}
void Finley_ElementFile_setElementDistribution(Finley_ElementFile* in, dim_t* distribution) {
  dim_t local_num_elements,e,num_elements=0, size;
  Paso_MPI_rank myRank;
  if (in == NULL) {
      distribution[0]=num_elements;
  } else {
      if (in->MPIInfo->size>1) {
         num_elements=0;
         myRank=in->MPIInfo->rank;
         size=in->MPIInfo->size;
         #pragma omp parallel private(local_num_elements)
         {
            local_num_elements=0;
            #pragma omp for private(e)
            for (e=0;e<in->numElements;e++) {
               if (in->Owner[e] == myRank) local_num_elements++;
            }
            #pragma omp critical
            num_elements+=local_num_elements;
         }
         #ifdef PASO_MPI
           MPI_Allgather(&num_elements,1,MPI_INT,distribution,1,MPI_INT,in->MPIInfo->comm);
         #else
           distribution[0]=num_elements;
         #endif
      } else {
        distribution[0]=in->numElements;
      }
  }
}

dim_t Finley_ElementFile_getGlobalNumElements(Finley_ElementFile* in) {
  dim_t size, *distribution=NULL, out, p;
  if (in == NULL) {
      return 0;
  } else {
    size=in->MPIInfo->size;
    distribution=TMPMEMALLOC(size,dim_t);
    Finley_ElementFile_setElementDistribution(in,distribution);
    out=0;
    for (p=0;p<size;++p) out+=distribution[p];
    TMPMEMFREE(distribution);
    return out;
  }
}
dim_t Finley_ElementFile_getMyNumElements(Finley_ElementFile* in) {
  dim_t size, *distribution=NULL, out;
  if (in == NULL) {
      return 0;
  } else {
    size=in->MPIInfo->size;
    distribution=TMPMEMALLOC(size,dim_t);
    Finley_ElementFile_setElementDistribution(in,distribution);
    out=distribution[in->MPIInfo->rank];
    TMPMEMFREE(distribution);
    return out;
  }

}
index_t Finley_ElementFile_getFirstElement(Finley_ElementFile* in){
  dim_t size, *distribution=NULL, out, p;
  if (in == NULL) {
      return 0;
  } else {
    size=in->MPIInfo->size;
    distribution=TMPMEMALLOC(size,dim_t);
    Finley_ElementFile_setElementDistribution(in,distribution);
    out=0;
    for (p=0;p<in->MPIInfo->rank;++p) out+=distribution[p];
    TMPMEMFREE(distribution);
    return out;
  }
}
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	/* Finley: ElementFile */
18
19	/* allocates an element file to hold elements of type id and with integration order order. */
20	/* use Finley_Mesh_allocElementTable to allocate the element table (Id,Nodes,Tag,Owner). */
21
22	/**************************************************************/
23
24	#include "ElementFile.h"
25
26	/**************************************************************/
27
28	Finley_ElementFile* Finley_ElementFile_alloc(ElementTypeId id, index_t order, index_t reduced_order, Paso_MPIInfo *MPIInfo)
29	{
30	extern Finley_RefElementInfo Finley_RefElement_InfoList[];
31	dim_t NQ, reduced_NQ;
32	Finley_ElementFile *out;
33
34	/* get the number of quadrature nodes needed to achieve integration order order: */
35
36	if (order<0) order=MAX(2*Finley_RefElement_InfoList[id].numOrder,0);
37	if (reduced_order<0) reduced_order=MAX(2*(Finley_RefElement_InfoList[id].numOrder-1),0);
38	NQ= Finley_RefElement_InfoList[id].getNumQuadNodes(order);
39	reduced_NQ= Finley_RefElement_InfoList[id].getNumQuadNodes(reduced_order);
40	if (! Finley_noError()) return NULL;
41
42	/* allocate the return value */
43
44	out=MEMALLOC(1,Finley_ElementFile);
45	if (Finley_checkPtr(out)) return NULL;
46	out->order = order;
47	out->reduced_order = reduced_order;
48	out->ReferenceElement=NULL;
49	out->LinearReferenceElement=NULL;
50	out->ReferenceElementReducedOrder=NULL;
51	out->LinearReferenceElementReducedOrder=NULL;
52	out->numElements=0;
53	out->Id=NULL;
54	out->Nodes=NULL;
55	out->Tag=NULL;
56	out->Color=NULL;
57	out->minColor=0;
58	out->maxColor=-1;
59	out->jacobeans=NULL;
60	out->jacobeans_reducedQ=NULL;
61	out->jacobeans_reducedS=NULL;
62	out->jacobeans_reducedS_reducedQ=NULL;
63
64	out->Owner=NULL;
65	out->numTagsInUse=0;
66	out->tagsInUse=NULL;
67
68	out->MPIInfo = Paso_MPIInfo_getReference( MPIInfo );
69
70	/* allocate the reference element: */
71
72	out->ReferenceElement=Finley_RefElement_alloc(id,NQ);
73	out->jacobeans=Finley_ElementFile_Jacobeans_alloc(out->ReferenceElement);
74	out->ReferenceElementReducedOrder=Finley_RefElement_alloc(id,reduced_NQ);
75	out->jacobeans_reducedQ=Finley_ElementFile_Jacobeans_alloc(out->ReferenceElementReducedOrder);
76
77	out->LinearReferenceElement=Finley_RefElement_alloc(Finley_RefElement_InfoList[id].LinearTypeId,NQ);
78	out->jacobeans_reducedS=Finley_ElementFile_Jacobeans_alloc(out->LinearReferenceElement);
79	out->LinearReferenceElementReducedOrder=Finley_RefElement_alloc(Finley_RefElement_InfoList[id].LinearTypeId,reduced_NQ);
80	out->jacobeans_reducedS_reducedQ=Finley_ElementFile_Jacobeans_alloc(out->LinearReferenceElementReducedOrder);
81
82	out->numNodes=out->ReferenceElement->Type->numNodes;
83
84	if (! Finley_noError()) {
85	Finley_ElementFile_free(out);
86	return NULL;
87	}
88	return out;
89	}
90
91	/* deallocates an element file: */
92
93	void Finley_ElementFile_free(Finley_ElementFile* in) {
94	if (in!=NULL) {
95	#ifdef Finley_TRACE
96	if (in->ReferenceElement!=NULL) printf("element file for %s is deallocated.\n",in->ReferenceElement->Type->Name);
97	#endif
98	Finley_ElementFile_freeTable(in);
99	Finley_RefElement_dealloc(in->ReferenceElement);
100	Finley_RefElement_dealloc(in->ReferenceElementReducedOrder);
101	Finley_RefElement_dealloc(in->LinearReferenceElement);
102	Finley_RefElement_dealloc(in->LinearReferenceElementReducedOrder);
103	Finley_ElementFile_Jacobeans_dealloc(in->jacobeans);
104	Finley_ElementFile_Jacobeans_dealloc(in->jacobeans_reducedS);
105	Finley_ElementFile_Jacobeans_dealloc(in->jacobeans_reducedQ);
106	Finley_ElementFile_Jacobeans_dealloc(in->jacobeans_reducedS_reducedQ);
107	Paso_MPIInfo_free( in->MPIInfo );
108	MEMFREE(in);
109	}
110	}
111	void Finley_ElementFile_setElementDistribution(Finley_ElementFile* in, dim_t* distribution) {
112	dim_t local_num_elements,e,num_elements=0, size;
113	Paso_MPI_rank myRank;
114	if (in == NULL) {
115	distribution[0]=num_elements;
116	} else {
117	if (in->MPIInfo->size>1) {
118	num_elements=0;
119	myRank=in->MPIInfo->rank;
120	size=in->MPIInfo->size;
121	#pragma omp parallel private(local_num_elements)
122	{
123	local_num_elements=0;
124	#pragma omp for private(e)
125	for (e=0;e<in->numElements;e++) {
126	if (in->Owner[e] == myRank) local_num_elements++;
127	}
128	#pragma omp critical
129	num_elements+=local_num_elements;
130	}
131	#ifdef PASO_MPI
132	MPI_Allgather(&num_elements,1,MPI_INT,distribution,1,MPI_INT,in->MPIInfo->comm);
133	#else
134	distribution[0]=num_elements;
135	#endif
136	} else {
137	distribution[0]=in->numElements;
138	}
139	}
140	}
141
142	dim_t Finley_ElementFile_getGlobalNumElements(Finley_ElementFile* in) {
143	dim_t size, *distribution=NULL, out, p;
144	if (in == NULL) {
145	return 0;
146	} else {
147	size=in->MPIInfo->size;
148	distribution=TMPMEMALLOC(size,dim_t);
149	Finley_ElementFile_setElementDistribution(in,distribution);
150	out=0;
151	for (p=0;p<size;++p) out+=distribution[p];
152	TMPMEMFREE(distribution);
153	return out;
154	}
155	}
156	dim_t Finley_ElementFile_getMyNumElements(Finley_ElementFile* in) {
157	dim_t size, *distribution=NULL, out;
158	if (in == NULL) {
159	return 0;
160	} else {
161	size=in->MPIInfo->size;
162	distribution=TMPMEMALLOC(size,dim_t);
163	Finley_ElementFile_setElementDistribution(in,distribution);
164	out=distribution[in->MPIInfo->rank];
165	TMPMEMFREE(distribution);
166	return out;
167	}
168
169	}
170	index_t Finley_ElementFile_getFirstElement(Finley_ElementFile* in){
171	dim_t size, *distribution=NULL, out, p;
172	if (in == NULL) {
173	return 0;
174	} else {
175	size=in->MPIInfo->size;
176	distribution=TMPMEMALLOC(size,dim_t);
177	Finley_ElementFile_setElementDistribution(in,distribution);
178	out=0;
179	for (p=0;p<in->MPIInfo->rank;++p) out+=distribution[p];
180	TMPMEMFREE(distribution);
181	return out;
182	}
183	}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision