finley/src/Assemble_LumpedSystem.c

/*
 ************************************************************
 *          Copyright 2006,2007 by ACcENULLNULL MNRF              *
 *                                                          *
 *              http://www.access.edu.au                    *
 *       Primary Business: Queensland, Australia            *
 *  Licensed under the Open NULLoftware License version 3.0    *
 *     http://www.opensource.org/licenses/osl-3.0.php       *
 *                                                          *
 ************************************************************
*/


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

/*    assembles the mass matrix in lumped form                */

/*    The coefficient D has to be defined on the integration points or not present. */

/*    lumpedMat has to be initialized before the routine is called. */

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

/*  Author: gross@access.edu.au */
/*  Version: $Id$ */

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

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


#define NEW_LUMPING /* */

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

void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D) 
{

  bool_t reducedIntegrationOrder=FALSE, expandedD;
  char error_msg[LenErrorMsg_MAX];
  Assemble_Parameters p;
  double time0;
  dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
  type_t funcspace;
  index_t color,*row_index=NULL;
  double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *D_p=NULL, *lumpedMat_p=NULL;
  register double rtmp, m_t, diagS;
  size_t len_EM_lumpedMat_size;

  Finley_resetError();

  if (nodes==NULL || elements==NULL) return;
  if (isEmpty(lumpedMat) || isEmpty(D)) return;
  if (isEmpty(lumpedMat) && !isEmpty(D)) { 
        Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficients are non-zero but no lumped matrix is given.");
        return;
  }
  funcspace=getFunctionSpaceType(D);
  /* check if all function spaces are the same */
  if (funcspace==FINLEY_ELEMENTS) {
       reducedIntegrationOrder=FALSE;
  } else if (funcspace==FINLEY_FACE_ELEMENTS)  {
       reducedIntegrationOrder=FALSE;
  } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
       reducedIntegrationOrder=TRUE;
  } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS)  {
       reducedIntegrationOrder=TRUE;
  } else {
       Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: assemblage failed because of illegal function space.");
  }
  if (! Finley_noError()) return;

  /* set all parameters in p*/
  Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
  if (! Finley_noError()) return;

  /* check if all function spaces are the same */

  if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
        sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
        Finley_setError(TYPE_ERROR,error_msg);
  }

  /*  check the dimensions: */
  
  if (p.numEqu==1) {
    if (!isEmpty(D)) {
       if (!isDataPointShapeEqual(D,0,dimensions)) {
          Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficient D, rank 0 expected.");
       }

    }
  } else {
    if (!isEmpty(D)) {
      dimensions[0]=p.numEqu;
      if (!isDataPointShapeEqual(D,1,dimensions)) {
          sprintf(error_msg,"Assemble_LumpedSystem: coefficient D, expected shape (%d,)",dimensions[0]);
          Finley_setError(TYPE_ERROR,error_msg);
      }
    }
  }

  if (Finley_noError()) {
    lumpedMat_p=getSampleData(lumpedMat,0);
    len_EM_lumpedMat=p.row_NN*p.numEqu;
    len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
    expandedD=isExpanded(D);
    S=p.row_jac->ReferenceElement->S;

    #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, m_t, diagS)
    {
       EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
       row_index=THREAD_MEMALLOC(p.row_NN,index_t);
       if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {
          if (p.numEqu == 1) {
             if (expandedD) {
                 #ifndef PASO_MPI
                 for (color=elements->minColor;color<=elements->maxColor;color++) {
                    /*  open loop over all elements: */
                    #pragma omp for private(e) schedule(static)
                    for(e=0;e<elements->numElements;e++){
                       if (elements->Color[e]==color) {
                 #else
                 {
                    for(e=0;e<elements->numElements;e++) {
                       {
                 #endif
                          Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
                          memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
                          D_p=getSampleData(D,e);
                          #ifdef NEW_LUMPING /* HRZ lumping */
                          /*
                           *           Number of PDEs: 1
                           *  D_p varies over element: True
                           */
                          m_t=0; /* mass of the element: m_t */
                          for (q=0;q<p.numQuad;q++) {
                              m_t+=Vol[q]*D_p[q];
                          }                           
                          diagS=0; /* diagonal sum: S */
                          for (s=0;s<p.row_NS;s++) {
                              rtmp=0;
                              for (q=0;q<p.numQuad;q++) {
                                  rtmp+=Vol[q]*D_p[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
                              }
                              EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
                              diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
                          }
                          /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
                          for (s=0;s<p.row_NS;s++) {
                              EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
                          }
                          #else /* row-sum lumping */
                          for (s=0;s<p.row_NS;s++) {
                              rtmp=0;
                              for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[q];
                              EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
                          }
                          #endif
                          for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
                          Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
                       } /* end color check */
                    } /* end element loop */
                  } /* end color loop */
             } else  {
                 #ifndef PASO_MPI
                 for (color=elements->minColor;color<=elements->maxColor;color++) {
                 /*  open loop over all elements: */
                 #pragma omp for private(e) schedule(static)
                 for(e=0;e<elements->numElements;e++){
                    if (elements->Color[e]==color) {
                 #else
                 {
                    for(e=0;e<elements->numElements;e++) {
                       {
                 #endif
                          Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
                          memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
                          D_p=getSampleData(D,e);
                          #ifdef NEW_LUMPING /* HRZ lumping */
                          /*
                           *           Number of PDEs: 1
                           *  D_p varies over element: False
                           */
                          m_t=0; /* mass of the element: m_t */
                          for (q=0;q<p.numQuad;q++) {
                              m_t+=Vol[q]*D_p[0];
                          }                           
                          diagS=0; /* diagonal sum: S */
                          for (s=0;s<p.row_NS;s++) {
                              rtmp=0;
                              for (q=0;q<p.numQuad;q++) {
                                  rtmp+=Vol[q]*D_p[0]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
                              }
                              EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
                              diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
                          }
                          /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
                          for (s=0;s<p.row_NS;s++) {
                              EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
                          }
                          #else /* row-sum lumping */
                          for (s=0;s<p.row_NS;s++) {
                              rtmp=0;
                              for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
                              EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
                          }
                          #endif
                          for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
                          Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
                      } /* end color check */
                    } /* end element loop */
                 } /* end color loop */
             }
          } else {
             if (expandedD) {
                 #ifndef PASO_MPI
                 for (color=elements->minColor;color<=elements->maxColor;color++) {
                    /*  open loop over all elements: */
                    #pragma omp for private(e) schedule(static)
                    for(e=0;e<elements->numElements;e++){
                       if (elements->Color[e]==color) {
                 #else
                 {
                    for(e=0;e<elements->numElements;e++) {
                       {
                 #endif
                          Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
                          memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
                          D_p=getSampleData(D,e);
                          #ifdef NEW_LUMPING /* HRZ lumping */
                          /*
                           *           Number of PDEs: Multiple
                           *  D_p varies over element: True
                           */
                          for (k=0;k<p.numEqu;k++) {
                              m_t=0; /* mass of the element: m_t */
                              for (q=0;q<p.numQuad;q++) {
                                  m_t+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)];
                              }                           
                              diagS=0; /* diagonal sum: S */
                              for (s=0;s<p.row_NS;s++) {
                                  rtmp=0;
                                  for (q=0;q<p.numQuad;q++) {
                                      rtmp+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
                                  }
                                  EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
                                  diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
                              }
                              /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
                              for (s=0;s<p.row_NS;s++) {
                                  EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
                              }
                          }
                          #else /* row-sum lumping */
                          for (s=0;s<p.row_NS;s++) {
                              for (k=0;k<p.numEqu;k++) {
                                  rtmp=0.;
                                  for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX2(k,q,p.numEqu)];
                                  EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
                              }
                          }
                          #endif
                          for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
                          Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
                       } /* end color check */
                    } /* end element loop */
                } /* end color loop */
             } else {
                 #ifndef PASO_MPI
                 for (color=elements->minColor;color<=elements->maxColor;color++) {
                    /*  open loop over all elements: */
                    #pragma omp for private(e) schedule(static)
                    for(e=0;e<elements->numElements;e++){
                       if (elements->Color[e]==color) {
                 #else
                 {
                    for(e=0;e<elements->numElements;e++) {
                       {
                 #endif
                          Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
                          memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
                          D_p=getSampleData(D,e);
                          #ifdef NEW_LUMPING /* HRZ lumping */
                          /*
                           *           Number of PDEs: Multiple
                           *  D_p varies over element: False
                           */
                          for (k=0;k<p.numEqu;k++) {
                              m_t=0; /* mass of the element: m_t */
                              for (q=0;q<p.numQuad;q++) {
                                  m_t+=Vol[q]*D_p[k];
                              }                           
                              diagS=0; /* diagonal sum: S */
                              for (s=0;s<p.row_NS;s++) {
                                  rtmp=0;
                                  for (q=0;q<p.numQuad;q++) {
                                      rtmp+=Vol[q]*D_p[k]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
                                  }
                                  EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
                                  diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
                              }
                              /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
                              for (s=0;s<p.row_NS;s++) {
                                  EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
                              }
                          }
                          #else /* row-sum lumping */
                          for (s=0;s<p.row_NS;s++) {
                              rtmp=0;
                              for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
                              for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
                          }
                          #endif
                          for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
                          Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
                       } /* end color check */
                    } /* end element loop */
                } /* end color loop */
             }
          }
       } /* end of pointer check */
       THREAD_MEMFREE(EM_lumpedMat);
       THREAD_MEMFREE(row_index);
    } /* end parallel region */
  }
  #ifdef Finley_TRACE
  printf("timing: assemblage lumped PDE: %.4e sec\n",Finley_timer()-time0);
  #endif
}
1	jgs	150	/*
2	elspeth	616	************************************************************
3	gross	1204	* Copyright 2006,2007 by ACcENULLNULL MNRF *
4	elspeth	616	* *
5			* http://www.access.edu.au *
6			* Primary Business: Queensland, Australia *
7	gross	1204	* Licensed under the Open NULLoftware License version 3.0 *
8	elspeth	616	* http://www.opensource.org/licenses/osl-3.0.php *
9			* *
10			************************************************************
11	jgs	150	*/
12	jgs	82
13	jgs	150
14	jgs	82	/**************************************************************/
15
16	gross	1204	/* assembles the mass matrix in lumped form */
17	jgs	82
18	gross	1204	/* The coefficient D has to be defined on the integration points or not present. */
19	jgs	82
20	gross	1204	/* lumpedMat has to be initialized before the routine is called. */
21	jgs	82
22			/**************************************************************/
23
24	jgs	150	/* Author: gross@access.edu.au */
25	btully	1220	/* Version: $Id$ */
26	jgs	82
27			/**************************************************************/
28
29			#include "Assemble.h"
30			#include "Util.h"
31			#ifdef _OPENMP
32			#include <omp.h>
33			#endif
34
35
36	btully	1220	#define NEW_LUMPING /* */
37
38	jgs	82	/**************************************************************/
39
40	gross	1204	void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D)
41			{
42	jgs	82
43	gross	1204	bool_t reducedIntegrationOrder=FALSE, expandedD;
44	jgs	150	char error_msg[LenErrorMsg_MAX];
45	gross	798	Assemble_Parameters p;
46	jgs	82	double time0;
47	gross	1204	dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
48	gross	1028	type_t funcspace;
49	gross	1204	index_t color,*row_index=NULL;
50			double S=NULL, EM_lumpedMat=NULL, Vol=NULL, D_p=NULL, *lumpedMat_p=NULL;
51	btully	1220	register double rtmp, m_t, diagS;
52	gross	1204	size_t len_EM_lumpedMat_size;
53	gross	798
54	jgs	150	Finley_resetError();
55	jgs	82
56			if (nodes==NULL \|\| elements==NULL) return;
57	gross	1204	if (isEmpty(lumpedMat) \|\| isEmpty(D)) return;
58			if (isEmpty(lumpedMat) && !isEmpty(D)) {
59			Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficients are non-zero but no lumped matrix is given.");
60			return;
61	gross	798	}
62	gross	1204	funcspace=getFunctionSpaceType(D);
63	jgs	82	/* check if all function spaces are the same */
64	gross	798	if (funcspace==FINLEY_ELEMENTS) {
65			reducedIntegrationOrder=FALSE;
66			} else if (funcspace==FINLEY_FACE_ELEMENTS) {
67			reducedIntegrationOrder=FALSE;
68	gross	1072	} else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
69			reducedIntegrationOrder=TRUE;
70			} else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
71			reducedIntegrationOrder=TRUE;
72	gross	798	} else {
73	gross	1204	Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: assemblage failed because of illegal function space.");
74	gross	798	}
75			if (! Finley_noError()) return;
76	jgs	82
77	gross	798	/* set all parameters in p*/
78	gross	1204	Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
79	gross	798	if (! Finley_noError()) return;
80
81			/* check if all function spaces are the same */
82
83	jgs	82	if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
84	gross	1204	sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
85	jgs	150	Finley_setError(TYPE_ERROR,error_msg);
86	jgs	82	}
87
88			/* check the dimensions: */
89
90	gross	1204	if (p.numEqu==1) {
91	jgs	82	if (!isEmpty(D)) {
92			if (!isDataPointShapeEqual(D,0,dimensions)) {
93	gross	1204	Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficient D, rank 0 expected.");
94	jgs	82	}
95	gross	1204
96	jgs	82	}
97			} else {
98			if (!isEmpty(D)) {
99			dimensions[0]=p.numEqu;
100	gross	1204	if (!isDataPointShapeEqual(D,1,dimensions)) {
101			sprintf(error_msg,"Assemble_LumpedSystem: coefficient D, expected shape (%d,)",dimensions[0]);
102	jgs	150	Finley_setError(TYPE_ERROR,error_msg);
103	jgs	82	}
104			}
105			}
106	gross	1204
107	jgs	150	if (Finley_noError()) {
108	gross	1204	lumpedMat_p=getSampleData(lumpedMat,0);
109			len_EM_lumpedMat=p.row_NN*p.numEqu;
110			len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
111			expandedD=isExpanded(D);
112			S=p.row_jac->ReferenceElement->S;
113
114	btully	1220	#pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, m_t, diagS)
115	gross	1204	{
116			EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
117			row_index=THREAD_MEMALLOC(p.row_NN,index_t);
118			if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {
119			if (p.numEqu == 1) {
120			if (expandedD) {
121			#ifndef PASO_MPI
122			for (color=elements->minColor;color<=elements->maxColor;color++) {
123			/* open loop over all elements: */
124			#pragma omp for private(e) schedule(static)
125			for(e=0;e<elements->numElements;e++){
126			if (elements->Color[e]==color) {
127			#else
128			{
129			for(e=0;e<elements->numElements;e++) {
130			{
131			#endif
132			Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
133			memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
134			D_p=getSampleData(D,e);
135	btully	1220	#ifdef NEW_LUMPING /* HRZ lumping */
136			/*
137			* Number of PDEs: 1
138			* D_p varies over element: True
139			*/
140			m_t=0; /* mass of the element: m_t */
141			for (q=0;q<p.numQuad;q++) {
142			m_t+=Vol[q]*D_p[q];
143			}
144			diagS=0; /* diagonal sum: S */
145	gross	1204	for (s=0;s<p.row_NS;s++) {
146			rtmp=0;
147	btully	1220	for (q=0;q<p.numQuad;q++) {
148			rtmp+=Vol[q]D_p[q]S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
149			}
150			EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
151			diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
152			}
153			/* rescale diagonals by m_t/diagS to ensure consistent mass over element */
154			for (s=0;s<p.row_NS;s++) {
155			EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
156			}
157			#else /* row-sum lumping */
158			for (s=0;s<p.row_NS;s++) {
159			rtmp=0;
160	gross	1204	for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]S[INDEX2(s,q,p.row_NS)]D_p[q];
161			EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
162			}
163	btully	1220	#endif
164	gross	1204	for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
165			Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
166			} /* end color check */
167			} /* end element loop */
168			} /* end color loop */
169			} else {
170			#ifndef PASO_MPI
171			for (color=elements->minColor;color<=elements->maxColor;color++) {
172			/* open loop over all elements: */
173			#pragma omp for private(e) schedule(static)
174			for(e=0;e<elements->numElements;e++){
175			if (elements->Color[e]==color) {
176			#else
177			{
178			for(e=0;e<elements->numElements;e++) {
179			{
180			#endif
181	btully	1220	Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
182			memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
183			D_p=getSampleData(D,e);
184			#ifdef NEW_LUMPING /* HRZ lumping */
185			/*
186			* Number of PDEs: 1
187			* D_p varies over element: False
188			*/
189			m_t=0; /* mass of the element: m_t */
190			for (q=0;q<p.numQuad;q++) {
191			m_t+=Vol[q]*D_p[0];
192			}
193			diagS=0; /* diagonal sum: S */
194			for (s=0;s<p.row_NS;s++) {
195			rtmp=0;
196			for (q=0;q<p.numQuad;q++) {
197			rtmp+=Vol[q]D_p[0]S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
198			}
199			EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
200			diagS+=EM_lumpedMat[INDEX2(0,s,p.numEqu)];
201			}
202			/* rescale diagonals by m_t/diagS to ensure consistent mass over element */
203			for (s=0;s<p.row_NS;s++) {
204			EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=m_t/diagS;
205			}
206			#else /* row-sum lumping */
207			for (s=0;s<p.row_NS;s++) {
208			rtmp=0;
209			for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
210			EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
211			}
212			#endif
213			for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
214			Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
215	gross	1204	} /* end color check */
216			} /* end element loop */
217			} /* end color loop */
218			}
219	gross	798	} else {
220	gross	1204	if (expandedD) {
221			#ifndef PASO_MPI
222			for (color=elements->minColor;color<=elements->maxColor;color++) {
223			/* open loop over all elements: */
224			#pragma omp for private(e) schedule(static)
225			for(e=0;e<elements->numElements;e++){
226			if (elements->Color[e]==color) {
227			#else
228			{
229			for(e=0;e<elements->numElements;e++) {
230			{
231			#endif
232			Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
233			memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
234			D_p=getSampleData(D,e);
235	btully	1220	#ifdef NEW_LUMPING /* HRZ lumping */
236			/*
237			* Number of PDEs: Multiple
238			* D_p varies over element: True
239			*/
240			for (k=0;k<p.numEqu;k++) {
241			m_t=0; /* mass of the element: m_t */
242			for (q=0;q<p.numQuad;q++) {
243			m_t+=Vol[q]*D_p[INDEX2(k,q,p.numEqu)];
244			}
245			diagS=0; /* diagonal sum: S */
246			for (s=0;s<p.row_NS;s++) {
247			rtmp=0;
248			for (q=0;q<p.numQuad;q++) {
249			rtmp+=Vol[q]D_p[INDEX2(k,q,p.numEqu)]S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
250			}
251			EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
252			diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
253			}
254			/* rescale diagonals by m_t/diagS to ensure consistent mass over element */
255			for (s=0;s<p.row_NS;s++) {
256			EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
257			}
258			}
259			#else /* row-sum lumping */
260	gross	1204	for (s=0;s<p.row_NS;s++) {
261			for (k=0;k<p.numEqu;k++) {
262			rtmp=0.;
263			for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]S[INDEX2(s,q,p.row_NS)]D_p[INDEX2(k,q,p.numEqu)];
264			EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
265			}
266			}
267	btully	1220	#endif
268	gross	1204	for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
269			Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
270			} /* end color check */
271			} /* end element loop */
272			} /* end color loop */
273			} else {
274			#ifndef PASO_MPI
275			for (color=elements->minColor;color<=elements->maxColor;color++) {
276			/* open loop over all elements: */
277			#pragma omp for private(e) schedule(static)
278			for(e=0;e<elements->numElements;e++){
279			if (elements->Color[e]==color) {
280			#else
281			{
282			for(e=0;e<elements->numElements;e++) {
283			{
284			#endif
285			Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
286			memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
287			D_p=getSampleData(D,e);
288	btully	1220	#ifdef NEW_LUMPING /* HRZ lumping */
289			/*
290			* Number of PDEs: Multiple
291			* D_p varies over element: False
292			*/
293			for (k=0;k<p.numEqu;k++) {
294			m_t=0; /* mass of the element: m_t */
295			for (q=0;q<p.numQuad;q++) {
296			m_t+=Vol[q]*D_p[k];
297			}
298			diagS=0; /* diagonal sum: S */
299			for (s=0;s<p.row_NS;s++) {
300			rtmp=0;
301			for (q=0;q<p.numQuad;q++) {
302			rtmp+=Vol[q]D_p[k]S[INDEX2(s,q,p.row_NS)]*S[INDEX2(s,q,p.row_NS)];
303			}
304			EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
305			diagS+=EM_lumpedMat[INDEX2(k,s,p.numEqu)];
306			}
307			/* rescale diagonals by m_t/diagS to ensure consistent mass over element */
308			for (s=0;s<p.row_NS;s++) {
309			EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=m_t/diagS;
310			}
311			}
312			#else /* row-sum lumping */
313	gross	1204	for (s=0;s<p.row_NS;s++) {
314			rtmp=0;
315			for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
316			for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
317			}
318	btully	1220	#endif
319	gross	1204	for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
320			Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
321			} /* end color check */
322			} /* end element loop */
323			} /* end color loop */
324			}
325	gross	798	}
326	gross	1204	} /* end of pointer check */
327			THREAD_MEMFREE(EM_lumpedMat);
328			THREAD_MEMFREE(row_index);
329			} /* end parallel region */
330	jgs	82	}
331	gross	1204	#ifdef Finley_TRACE
332			printf("timing: assemblage lumped PDE: %.4e sec\n",Finley_timer()-time0);
333			#endif
334	jgs	82	}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision