Small changes for assembling.

AMDiS/libtool

@@ -82,13 +82,13 @@ AR="ar"
 AR_FLAGS="cru"
 
 # A C compiler.
-LTCC="/usr/lib/openmpi/1.2.7-gcc//bin/mpicc"
+LTCC="gcc"
 
 # LTCC compiler flags.
 LTCFLAGS="-g -O2"
 
 # A language-specific compiler.
-CC="/usr/lib/openmpi/1.2.7-gcc//bin/mpicc"
+CC="gcc"
 
 # Is the compiler the GNU C compiler?
 with_gcc=yes
@@ -174,7 +174,7 @@ dlopen_self=unknown
 dlopen_self_static=unknown
 
 # Compiler flag to prevent dynamic linking.
-link_static_flag=""
+link_static_flag="-static"
 
 # Compiler flag to turn off builtin functions.
 no_builtin_flag=" -fno-builtin"
@@ -6801,13 +6801,13 @@ AR="ar"
 AR_FLAGS="cru"
 
 # A C compiler.
-LTCC="/usr/lib/openmpi/1.2.7-gcc//bin/mpicc"
+LTCC="gcc"
 
 # LTCC compiler flags.
 LTCFLAGS="-g -O2"
 
 # A language-specific compiler.
-CC="/usr/lib/openmpi/1.2.7-gcc//bin/mpiCC"
+CC="g++"
 
 # Is the compiler the GNU C compiler?
 with_gcc=yes
@@ -6893,7 +6893,7 @@ dlopen_self=unknown
 dlopen_self_static=unknown
 
 # Compiler flag to prevent dynamic linking.
-link_static_flag=""
+link_static_flag="-static"
 
 # Compiler flag to turn off builtin functions.
 no_builtin_flag=" -fno-builtin"
@@ -6960,11 +6960,11 @@ predeps=""
 
 # Dependencies to place after the objects being linked to create a
 # shared library.
-postdeps="-lmpi_cxx -lmpi -lopen-rte -lopen-pal -ldl -lnsl -lutil -ldl -lstdc++ -lm -lgcc_s -lpthread -lc -lgcc_s"
+postdeps="-lstdc++ -lm -lgcc_s -lc -lgcc_s"
 
 # The library search path used internally by the compiler when linking
 # a shared library.
-compiler_lib_search_path=`echo "-L/usr/lib/openmpi/1.2.7-gcc/lib -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2/../../.." | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+compiler_lib_search_path=`echo "-L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2/../../.." | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
 
 # Method to check whether dependent libraries are shared objects.
 deplibs_check_method="pass_all"
@@ -7109,7 +7109,7 @@ AR="ar"
 AR_FLAGS="cru"
 
 # A C compiler.
-LTCC="/usr/lib/openmpi/1.2.7-gcc//bin/mpicc"
+LTCC="gcc"
 
 # LTCC compiler flags.
 LTCFLAGS="-g -O2"

AMDiS/src/FirstOrderAssembler.cc

@@ -199,10 +199,12 @@ namespace AMDiS {
 
       const double *phi = phiFast->getPhi(iq);
       grdPsi = psiFast->getGradient(iq);
+      double factor = quadrature->getWeight(iq);
 
-      for (int i = 0; i < nRow; i++)
+      for (int i = 0; i < nRow; i++) {
 	for (int j = 0; j < nCol; j++)
-	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j];
+	  mat[i][j] += factor * (Lb[iq] * (*grdPsi)[i]) * phi[j];
+      }
     }
   }
 
@@ -358,6 +360,7 @@ namespace AMDiS {
 
     for (int iq = 0; iq < nPoints; iq++)
       Lb[iq].set(0.0);
+
     for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
       (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
   
@@ -367,9 +370,11 @@ namespace AMDiS {
       const double *psi = psiFast->getPhi(iq);
       grdPhi = phiFast->getGradient(iq);
 
-      for (int i = 0; i < nRow; i++)
+      for (int i = 0; i < nRow; i++) {
+	double factor = quadrature->getWeight(iq) * psi[i];
 	for (int j = 0; j < nCol; j++)
-	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPhi)[j]) * psi[i];
+	  mat[i][j] += factor * (Lb[iq] * (*grdPhi)[j]);
+      }
     }
   }
 

AMDiS/src/Lagrange.cc

@@ -813,19 +813,23 @@ namespace AMDiS {
   {
     FUNCNAME("Lagrange::interpol()");
 
-    static double *inter = NULL;
-    static int inter_size = 0;
+    static double* localVec = NULL;
+    static int localVecSize = 0;
 
-    double *rvec =  NULL;
+    double *rvec = NULL;
 
     if (vec) {
       rvec = vec;
     } else {
-      if (inter) 
-	delete [] inter;
-      inter = new double[nBasFcts];
-      inter_size = nBasFcts;
-      rvec = inter;
+      if (localVec && nBasFcts > localVecSize)  {
+	delete [] localVec;
+	localVec = new double[nBasFcts]; 
+      }
+      if (!localVec)
+	localVec = new double[nBasFcts]; 
+
+      localVecSize = nBasFcts;
+      rvec = localVec;
     } 
 
     WorldVector<double> x;
@@ -1000,57 +1004,46 @@ namespace AMDiS {
     ERROR_EXIT("not yet\n");
   }
 
-  void Lagrange::l2ScpFctBas(Quadrature *q,
+  void Lagrange::l2ScpFctBas(Quadrature *quad,
 			     AbstractFunction<double, WorldVector<double> >* f,
 			     DOFVector<double>* fh)
   {
     FUNCNAME("Lagrange::l2ScpFctBas()");
 
-    TraverseStack stack;
-    double *wdetf_qp = NULL;
-    double val, det;
-    WorldVector<double> x;
-    int iq, j;
-    const DegreeOfFreedom  *dof;
-    const Parametric *parametric;
-
     TEST_EXIT_DBG(fh)("no DOF_REAL_VEC fh\n");
     TEST_EXIT_DBG(fh->getFESpace())
       ("no fe_space in DOF_REAL_VEC %s\n", fh->getName().c_str());
     TEST_EXIT_DBG(fh->getFESpace()->getBasisFcts() == this)
       ("wrong basis fcts for fh\n");
+    TEST_EXIT_DBG(!fh->getFESpace()->getMesh()->getParametric())
+      ("Not yet implemented!");
 
-    Mesh* mesh = fh->getFESpace()->getMesh();
-    int n_phi = nBasFcts;
-
-    if (!q)
-      q = Quadrature::provideQuadrature(dim, 2 * degree - 2);
+    if (!quad)
+      quad = Quadrature::provideQuadrature(dim, 2 * degree - 2);
 
     const FastQuadrature *quad_fast = 
-      FastQuadrature::provideFastQuadrature(this, *q, INIT_PHI);
-
-    if ((parametric = mesh->getParametric()))
-      ERROR_EXIT("not yet implemented\n");
-    else
-      wdetf_qp = new double[q->getNumPoints()];
+      FastQuadrature::provideFastQuadrature(this, *quad, INIT_PHI);
+    double *wdetf_qp = new double[quad->getNumPoints()];
+    int nPoints = quad->getNumPoints();
+    DOFAdmin *admin = fh->getFESpace()->getAdmin();
+    WorldVector<double> x;
 
-    ElInfo *el_info = stack.traverseFirst(mesh, -1, 
+    TraverseStack stack;
+    ElInfo *el_info = stack.traverseFirst(fh->getFESpace()->getMesh(), -1, 
 					  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
-
-    int numPoints = q->getNumPoints();
-     
     while (el_info) {
-      dof = getLocalIndices(el_info->getElement(), (fh->getFESpace()->getAdmin()), NULL);      
-      det = el_info->getDet();
+      const DegreeOfFreedom  *dof = getLocalIndices(el_info->getElement(), admin, NULL);
+      double det = el_info->getDet();
 
-      for (iq = 0; iq < numPoints; iq++) {
-	el_info->coordToWorld(q->getLambda(iq), x);
-	wdetf_qp[iq] = q->getWeight(iq)*det*((*f)(x));
+      for (int iq = 0; iq < nPoints; iq++) {
+	el_info->coordToWorld(quad->getLambda(iq), x);
+	wdetf_qp[iq] = quad->getWeight(iq) * det * ((*f)(x));
       }
       
-      for (j = 0; j < n_phi; j++) {
-	for (val = iq = 0; iq < numPoints; iq++)
-	  val += quad_fast->getPhi(iq,j)*wdetf_qp[iq];
+      for (int j = 0; j < nBasFcts; j++) {
+	double val = 0.0;
+	for (int iq = 0; iq < nPoints; iq++)
+	  val += quad_fast->getPhi(iq, j) * wdetf_qp[iq];
 
 	(*fh)[dof[j]] += val;
       }
@@ -1096,9 +1089,12 @@ namespace AMDiS {
     int dim = drv->getFESpace()->getMesh()->getDim();
     int n0 = drv->getFESpace()->getAdmin()->getNumberOfPreDOFs(VERTEX);
     Element *el = list->getElement(0);
-    DegreeOfFreedom dof0 = el->getDOF(0, n0);           /* 1st endpoint of refinement edge */
-    DegreeOfFreedom dof1 = el->getDOF(1, n0);           /* 2nd endpoint of refinement edge */
-    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(dim, n0);  /*      newest vertex is DIM */
+    // 1st endpoint of refinement edge
+    DegreeOfFreedom dof0 = el->getDOF(0, n0);
+    // 2nd endpoint of refinement edge          
+    DegreeOfFreedom dof1 = el->getDOF(1, n0);
+    // newest vertex is DIM 
+    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(dim, n0); 
     (*drv)[dof_new] = 0.5 * ((*drv)[dof0] + (*drv)[dof1]);
   }
 
@@ -1121,7 +1117,8 @@ namespace AMDiS {
     /*  newest vertex of child[0] and child[1]                                  */
     /****************************************************************************/
 
-    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node+1, n0);  /*      newest vertex is DIM */
+    // newest vertex is DIM
+    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node + 1, n0);
     (*drv)[cdof] = (*drv)[pdof[2]];
 
     node = drv->getFESpace()->getMesh()->getNode(CENTER);        

AMDiS/src/Operator.cc

@@ -245,22 +245,6 @@ namespace AMDiS {
     }      
   }
 
-  void OperatorTerm::lb(const DimVec<WorldVector<double> >& Lambda,
-			const WorldVector<double>& b,
-			DimVec<double>& Lb,
-			double factor)
-  {
-    int dim = Lb.getSize() - 1;
-    const int dimOfWorld = Global::getGeo(WORLD);
-
-    for (int i = 0; i <= dim; i++) {
-      double val = 0.0;
-      for (int j = 0; j < dimOfWorld; j++)
-	val += Lambda[i][j] * b[j];
-      Lb[i] += val * factor;
-    }    
-  }
-
   Operator::Operator(Flag operatorType,
 		     const FiniteElemSpace *row,
 		     const FiniteElemSpace *col)

AMDiS/src/Operator.h

@@ -166,27 +166,40 @@ namespace AMDiS {
 		     double factor);
 
     /// Evaluation of \f$ \Lambda \cdot b\f$.
-    static void lb(const DimVec<WorldVector<double> >& Lambda,
-		   const WorldVector<double>& b,
-		   DimVec<double>& Lb,
-		   double factor);
-
-    /// Evaluation of \f$ \Lambda \cdot b\f$ if b contains the value 1.0 in each component.
-    static void l1(const DimVec<WorldVector<double> >& Lambda,
-		   DimVec<double>& Lb,
-		   double factor)
+    static inline void lb(const DimVec<WorldVector<double> >& Lambda,
+			  const WorldVector<double>& b,
+			  DimVec<double>& Lb,
+			  double factor)
+    {
+      const int dimOfWorld = Global::getGeo(WORLD);
+      
+      for (int i = 0; i <= dimOfWorld; i++) {
+	double val = 0.0;
+	
+	for (int j = 0; j < dimOfWorld; j++)
+	  val += Lambda[i][j] * b[j];
+	
+	Lb[i] += val * factor;
+      }    
+    }
+
+    /** \brief
+     * Evaluation of \f$ \Lambda \cdot b\f$ if b contains the value 1.0 in 
+     * each component.
+     */
+    static inline void l1(const DimVec<WorldVector<double> >& Lambda,
+			  DimVec<double>& Lb,
+			  double factor)
     {
-      int dim = Lb.getSize() - 1;
       const int dimOfWorld = Global::getGeo(WORLD);
 
-      for (int i = 0; i <= dim; i++) {
+      for (int i = 0; i <= dimOfWorld; i++) {
 	double val = 0.0;
       
 	for (int j = 0; j < dimOfWorld; j++)
 	  val += Lambda[i][j];
 
-	val *= factor;
-	Lb[i] += val;
+	Lb[i] += val * factor;
       }    
     }