DENDRO: /net/www/csela/dendro/releases/Dendro/examples/src/backend/matProp.C Source File

00001 
00007 #include "omg.h"
00008 #include "oda.h" 
00009 #include "omgJac.h"
00010 #include "parUtils.h"
00011 
00012 #define SQUARE(x) ((x)*(x))
00013 
00014 #define ASSIGN_MAT_PROP_FROM_PTS_ELEM_BLOCK {\
00015   unsigned int idx = da->curr();\
00016   unsigned int lev = da->getLevel(idx);\
00017   Point pt = da->getCurrentOffset();\
00018   ot::TreeNode me(pt.xint(),pt.yint(),pt.zint(),lev,3,maxD);\
00019   matPropArr[2*idx] = lapBase; /*Laplacian's coefficient*/\
00020   matPropArr[2*idx+1] = massBase; /*Mass Matrix's coefficient*/\
00021   if(ptsCtr < (lapJump.size())) {\
00022     ot::TreeNode lastVisitedPt;\
00023     /*Assuming pts are sorted in Morton ordering*/\
00024     /*Assuming pts and local octants are aligned*/\
00025     bool repeat = false;\
00026     do {\
00027       unsigned int ptX = static_cast<unsigned int>(\
00028           pts[(3*ptsCtr)]*static_cast<double>(1u << (maxD-1)));\
00029       unsigned int ptY = static_cast<unsigned int>(\
00030           pts[(3*ptsCtr)+1]*static_cast<double>(1u << (maxD-1)));\
00031       unsigned int ptZ = static_cast<unsigned int>(\
00032           pts[(3*ptsCtr)+2]*static_cast<double>(1u << (maxD-1)));\
00033       ot::TreeNode currPt(ptX,ptY,ptZ,maxD,3,maxD);\
00034       repeat = false;\
00035       if( (currPt < me) && (ptsCtr < (lapJump.size()-1)) ) {\
00036         repeat = true;\
00037         ptsCtr++;\
00038       } else {\
00039         lastVisitedPt = currPt;\
00040       }\
00041     } while(repeat);\
00042     /*There could be more than 1 pt. per octant*/\
00043     while( (ptsCtr < lapJump.size()) && (me.isAncestor(lastVisitedPt) ||\
00044           (me == lastVisitedPt))) {\
00045       matPropArr[2*idx] += lapJump[ptsCtr];\
00046       ptsCtr++;\
00047       if(ptsCtr < lapJump.size()) {\
00048         unsigned int ptX = static_cast<unsigned int>(\
00049             pts[(3*ptsCtr)]*static_cast<double>(1u << (maxD-1)));\
00050         unsigned int ptY = static_cast<unsigned int>(\
00051             pts[(3*ptsCtr)+1]*static_cast<double>(1u << (maxD-1)));\
00052         unsigned int ptZ = static_cast<unsigned int>(\
00053             pts[(3*ptsCtr)+2]*static_cast<double>(1u << (maxD-1)));\
00054         ot::TreeNode currPt(ptX,ptY,ptZ,maxD,3,maxD);\
00055         lastVisitedPt = currPt;\
00056       }\
00057     }\
00058   }\
00059   if(matPropArr[2*idx] > maxCoeff) {\
00060     maxCoeff = matPropArr[2*idx];\
00061   }\
00062   if(matPropArr[2*idx] < minCoeff) {\
00063     minCoeff = matPropArr[2*idx];\
00064   }\
00065 }
00066 
00067 #define ASSIGN_MAT_PROP_1_CUBE_ELEM_BLOCK {\
00068   Point pt = da->getCurrentOffset();\
00069   unsigned int idx = da->curr();\
00070   unsigned int lev = da->getLevel(idx);\
00071   double h = ((double)(1u << (maxD - lev)))/((double)(1u << (maxD-1)));\
00072   double x = ((double)(pt.xint()))/((double)(1u << (maxD-1)));\
00073   double y = ((double)(pt.yint()))/((double)(1u << (maxD-1)));\
00074   double z = ((double)(pt.zint()))/((double)(1u << (maxD-1)));\
00075   matPropArr[2*idx] = 1.0; /*Laplacian's coefficient*/\
00076   matPropArr[2*idx+1] = massBase; /*Mass Matrix's coefficient*/\
00077   double minX1 = 0.5-h;\
00078   double minY1 = 0.0;\
00079   double minZ1 = 0.0;\
00080   double maxX1 = 0.5+h;\
00081   double maxY1 = 1.0;\
00082   double maxZ1 = 1.0;\
00083   double minX2 = x;\
00084   double minY2 = y;\
00085   double minZ2 = z;\
00086   double maxX2 = (x+h);\
00087   double maxY2 = (y+h);\
00088   double maxZ2 = (z+h);\
00089   bool oneInTwo = ( (minX1 >= minX2) && (minX1 < maxX2) &&  (minY1 >= minY2) && (minY1 < maxY2) &&  (minZ1 >= minZ2) && (minZ1 < maxZ2) );\
00090   bool twoInOne = ( (minX2 >= minX1) && (minX2 < maxX1) &&  (minY2 >= minY1) && (minY2 < maxY1) &&  (minZ2 >= minZ1) && (minZ2 < maxZ1) );\
00091   if( oneInTwo || twoInOne ) {\
00092     double maxOfMinsX = ( (minX2 > minX1) ? minX2 : minX1 );\
00093     double maxOfMinsY = ( (minY2 > minY1) ? minY2 : minY1 );\
00094     double maxOfMinsZ = ( (minZ2 > minZ1) ? minZ2 : minZ1 );\
00095     double minOfMaxsX = ( (maxX2 < maxX1) ? maxX2 : maxX1 );\
00096     double minOfMaxsY = ( (maxY2 < maxY1) ? maxY2 : maxY1 );\
00097     double minOfMaxsZ = ( (maxZ2 < maxZ1) ? maxZ2 : maxZ1 );\
00098     double intersectionVolume = ((minOfMaxsX - maxOfMinsX)*(minOfMaxsY - maxOfMinsY)*(minOfMaxsZ - maxOfMinsZ))/(h*h*h);\
00099     if(intersectionVolume <= 0.0) {\
00100       std::cout<<"intesection Volume: "<<intersectionVolume<<" x: "<<x<<" y: "<<y<<" z: "<<z<<" h: "<<h<<" oneInTwo: "<<oneInTwo<<" twoInOne: "<<twoInOne<<std::endl;\
00101       assert(false);\
00102     }\
00103     matPropArr[2*idx] += lapFac*intersectionVolume; /*Laplacian's coefficient*/\
00104     matPropArr[2*idx+1]+= massFac*intersectionVolume;/*Mass Matrix's coefficient*/\
00105   }\
00106   if(matPropArr[2*idx] > maxCoeff) {\
00107     maxCoeff = matPropArr[2*idx];\
00108   }\
00109   if(matPropArr[2*idx] < minCoeff) {\
00110     minCoeff = matPropArr[2*idx];\
00111   }\
00112 }
00113 
00114 #define ASSIGN_MAT_PROP_MULTIPLE_CUBES_ELEM_BLOCK {\
00115   unsigned int idx = da->curr();\
00116   matPropArr[2*idx] = 1.0; /*Laplacian's coefficient*/\
00117   matPropArr[2*idx+1] = massBase; /*Mass Matrix's coefficient*/\
00118   unsigned int procElementSize = da->getElementSize();\
00119   unsigned int numCubesInterval = procElementSize/numCubes;\
00120   if( (idx % numCubesInterval) == 0 ) {\
00121     matPropArr[2*idx] += lapFac; /*Laplacian's coefficient*/\
00122     matPropArr[2*idx+1] +=  massFac; /*Mass Matrix's coefficient*/\
00123   }\
00124   if(matPropArr[2*idx] > maxCoeff) {\
00125     maxCoeff = matPropArr[2*idx];\
00126   }\
00127   if(matPropArr[2*idx] < minCoeff) {\
00128     minCoeff = matPropArr[2*idx];\
00129   }\
00130 }
00131 
00132 #define ASSIGN_MAT_PROP_CHECKER_BOARD_ELEM_BLOCK {\
00133   Point pt = da->getCurrentOffset();\
00134   unsigned int idx = da->curr();\
00135   unsigned int lev = da->getLevel(idx);\
00136   double h = ((double)(1u << (maxD - lev)))/((double)(1u << (maxD-1)));\
00137   double x = ((double)(pt.xint()))/((double)(1u << (maxD-1)));\
00138   double y = ((double)(pt.yint()))/((double)(1u << (maxD-1)));\
00139   double z = ((double)(pt.zint()))/((double)(1u << (maxD-1)));\
00140   matPropArr[2*idx+1] = massBase; /*Mass Matrix's coefficient*/\
00141   if(z >= 0.0 && z < 0.5 ) {\
00142     if(x >= 0.0 && x < 0.5 && y >= 0.0 && y < 0.5) {\
00143       matPropArr[2*idx] = lapFac; /*Laplacian's coefficient*/\
00144     } else if(x >= 0.5 && x < 1.0 && y >= 0.0 && y < 0.5) {\
00145       matPropArr[2*idx] = 1.0; /*Laplacian's coefficient*/\
00146     } else if(x >= 0.0 && x < 0.5 && y >= 0.5 && y < 1.0) {\
00147       matPropArr[2*idx] = 1.0; /*Laplacian's coefficient*/\
00148     } else {\
00149       matPropArr[2*idx] = lapFac; /*Laplacian's coefficient*/\
00150     }\
00151   } else {\
00152     if(x >= 0.0 && x < 0.5 && y >= 0.0 && y < 0.5) {\
00153       matPropArr[2*idx] = 1.0; /*Laplacian's coefficient*/\
00154     } else if(x >= 0.5 && x < 1.0 && y >= 0.0 && y < 0.5) {\
00155       matPropArr[2*idx] = lapFac; /*Laplacian's coefficient*/\
00156     } else if(x >= 0.0 && x < 0.5 && y >= 0.5 && y < 1.0) {\
00157       matPropArr[2*idx] = lapFac; /*Laplacian's coefficient*/\
00158     } else {\
00159       matPropArr[2*idx] = 1.0; /*Laplacian's coefficient*/\
00160     }\
00161   }\
00162   if(matPropArr[2*idx] > maxCoeff) {\
00163     maxCoeff = matPropArr[2*idx];\
00164   }\
00165   if(matPropArr[2*idx] < minCoeff) {\
00166     minCoeff = matPropArr[2*idx];\
00167   }\
00168 }
00169 
00170 #define ASSIGN_MAT_PROP_ANALYTIC_FN_ELEM_BLOCK {\
00171   Point pt = da->getCurrentOffset();\
00172   unsigned int idx = da->curr();\
00173   unsigned int lev = da->getLevel(idx);\
00174   double h = ((double)(1u << (maxD - lev)))/((double)(1u << (maxD-1)));\
00175   double x = ((double)(pt.xint()))/((double)(1u << (maxD-1)));\
00176   double y = ((double)(pt.yint()))/((double)(1u << (maxD-1)));\
00177   double z = ((double)(pt.zint()))/((double)(1u << (maxD-1)));\
00178   matPropArr[2*idx+1] = massBase; /*Mass Matrix's coefficient*/\
00179   /*Laplacian's coefficient*/\
00180   matPropArr[2*idx] = 1.0 +\
00181   (lapFac*(SQUARE(cos(lapFreq*M_PI*(x+(0.5*h))))+\
00182            SQUARE(cos(lapFreq*M_PI*(y+(0.5*h))))+\
00183            SQUARE(cos(lapFreq*M_PI*(z+(0.5*h))))));\
00184   if(matPropArr[2*idx] > maxCoeff) {\
00185     maxCoeff = matPropArr[2*idx];\
00186   }\
00187   if(matPropArr[2*idx] < minCoeff) {\
00188     minCoeff = matPropArr[2*idx];\
00189   }\
00190 }
00191 
00192 #define ASSIGN_MAT_PROP_FINE_TO_COARSE_ELEM_BLOCK {\
00193   /*The fine loop is always Writable, but the coarse loop*/\
00194   /*could be Independent or W_Dependent. Hence the fine counter must*/\
00195   /*be incremented properly to align with the coarse.*/\
00196   unsigned int idxC= da->curr();\
00197   Point Cpt = da->getCurrentOffset();\
00198   assert(daf->curr() < daf->end<ot::DA_FLAGS::WRITABLE>());\
00199   while(daf->getCurrentOffset() != Cpt) {\
00200     daf->next<ot::DA_FLAGS::WRITABLE>();\
00201     assert(daf->curr() < daf->end<ot::DA_FLAGS::WRITABLE>());\
00202   }\
00203   if(daf->getLevel(daf->curr()) == da->getLevel(idxC)) {\
00204     /*The coarse and fine elements are the same,*/\
00205     assert(daf->curr() < daf->end<ot::DA_FLAGS::WRITABLE>());\
00206     unsigned int idxF = daf->curr();\
00207     matPropArr[2*idxC] = fMatArr[2*idxF];\
00208     matPropArr[2*idxC+1] = fMatArr[2*idxF+1];\
00209     daf->next<ot::DA_FLAGS::WRITABLE>();\
00210   }else {\
00211     double cLapVal = 0.0;\
00212     double cMassVal = 0.0;\
00213     for(unsigned char cNumFine = 0; cNumFine < 8; cNumFine++) {\
00214       /*The coarse and fine elements are NOT the same. */\
00215       /*Loop over each of the 8 children of the coarse element.*/\
00216       /*These are the underlying fine elements.*/\
00217       assert(daf->curr() < daf->end<ot::DA_FLAGS::WRITABLE>());\
00218       unsigned int idxF = daf->curr();\
00219       cLapVal += fMatArr[2*idxF];\
00220       cMassVal += fMatArr[2*idxF+1];\
00221       daf->next<ot::DA_FLAGS::WRITABLE>();\
00222     }\
00223     matPropArr[2*idxC] = (cLapVal/8.0);\
00224     matPropArr[2*idxC+1] = (cMassVal/8.0);\
00225   }\
00226   if(matPropArr[2*idxC] > maxCoeff) {\
00227     maxCoeff = matPropArr[2*idxC];\
00228   }\
00229   if(matPropArr[2*idxC] < minCoeff) {\
00230     minCoeff = matPropArr[2*idxC];\
00231   }\
00232 }
00233 
00234 #define CHK_AND_SCATTER_FINE_TO_COARSE_BLOCK {\
00235   if( damg_i->da_aux == NULL ) {\
00236     daf = da;\
00237     fMatPropVec = matPropVecPtr;\
00238     changedPartition  = false;\
00239   }else {\
00240     daf = damg_i->da_aux;\
00241     /*Need to Scatter Values*/\
00242     fMatPropVec = new std::vector<double>;\
00243     changedPartition  = true;\
00244     /*elemental - non-ghosted*/\
00245     std::vector<double> tmpVec1;\
00246     da->createVector<double>(tmpVec1,true,false,2);\
00247     double *vec1Arr = NULL;\
00248     /*Elemental,non-Ghosted,Write,2 Dof.*/\
00249     da->vecGetBuffer<double>(tmpVec1,vec1Arr,true,false,false,2);\
00250     matPropArr = NULL;\
00251     /*Elemental,Ghosted,Read-only,2 Dof.*/\
00252     da->vecGetBuffer<double>((*matPropVecPtr), matPropArr, true, true, true, 2);\
00253     if(da->iAmActive()) {\
00254       for(da->init<ot::DA_FLAGS::WRITABLE>();\
00255           da->curr() < da->end<ot::DA_FLAGS::WRITABLE>();\
00256           da->next<ot::DA_FLAGS::WRITABLE>()) {\
00257         unsigned int idx = da->curr();\
00258         vec1Arr[2*idx] = matPropArr[2*idx];\
00259         vec1Arr[2*idx+1] = matPropArr[2*idx+1];\
00260       }\
00261     }\
00262     da->vecRestoreBuffer<double>((*matPropVecPtr), matPropArr, true, true, true, 2);\
00263     da->vecRestoreBuffer<double>(tmpVec1,vec1Arr,true,false,false,2);\
00264     par::scatterValues<double>(tmpVec1, (*fMatPropVec), (2*(daf->getElementSize())), da->getComm());\
00265     tmpVec1.clear();\
00266   }\
00267 }
00268 
00269 #define SET_SINGLE_LEVEL_MAT_PROP_BLOCK(ELEM_MAT_PROP_BLOCK, EXTRA_INITS)  {\
00270   damg_i->user = ctx;\
00271   da = damg_i->da;\
00272   comm = da->getCommActive();\
00273   /*Elem,Ghosted, 2-dof vector.*/\
00274   /*Note: I am creating a ghosted vector only*/\
00275   /*because the mat-vec will need it.*/\
00276   /*So this way, I can avoid mallocs inside the mat-vec.*/\
00277   da->createVector<double>((*matPropVecPtr), true, true, 2);\
00278   for(unsigned int i = 0; i < matPropVecPtr->size(); i++) {\
00279     (*matPropVecPtr)[i] = 0.0;\
00280   }\
00281   /*Elemental,Ghosted,Write,2 Dof.*/\
00282   da->vecGetBuffer<double>((*matPropVecPtr), matPropArr,\
00283       true, true, false, 2);\
00284   maxD =  da->getMaxDepth();\
00285   if(da->iAmActive()) {\
00286     MPI_Comm_rank(comm,&rank);\
00287     double maxCoeff = 0.0;\
00288     double minCoeff = 1.0e+8;\
00289     double globalMaxCoeff;\
00290     double globalMinCoeff;\
00291     /*Any independent element has all its 8*/\
00292     /*indices pointing to elements the same processor owns.*/\
00293     /*So, even if any independent element is sent to*/\
00294     /*another processor as a ghost (This is rare, but might*/\
00295     /*happen due to a-priori communication and partial*/\
00296     /*second ring communication), it would have to*/\
00297     /*be FOREIGN on that processor. So overlapping comm*/\
00298     /*and comp is not a problem.*/\
00299     EXTRA_INITS\
00300     for(da->init<ot::DA_FLAGS::W_DEPENDENT>();\
00301         da->curr() < da->end<ot::DA_FLAGS::W_DEPENDENT>();\
00302         da->next<ot::DA_FLAGS::W_DEPENDENT>()) {\
00303       ELEM_MAT_PROP_BLOCK\
00304     }\
00305     da->ReadFromGhostElemsBegin<double>(matPropArr, 2);\
00306     EXTRA_INITS\
00307     for(da->init<ot::DA_FLAGS::INDEPENDENT>();\
00308         da->curr() < da->end<ot::DA_FLAGS::INDEPENDENT>();\
00309         da->next<ot::DA_FLAGS::INDEPENDENT>()) {\
00310       ELEM_MAT_PROP_BLOCK\
00311     }\
00312     da->ReadFromGhostElemsEnd<double>(matPropArr);\
00313     par::Mpi_Reduce<double>(&maxCoeff, &globalMaxCoeff, 1, MPI_MAX, 0, comm);\
00314     par::Mpi_Reduce<double>(&minCoeff, &globalMinCoeff, 1, MPI_MIN, 0, comm);\
00315     if(!rank) {\
00316       std::cout<<"Level: "<<(nlevels-(damg_i->nlevels))<<\
00317       " Max Lap. Coeff: "<<globalMaxCoeff<<\
00318       " Min Lap. Coeff: "<<globalMinCoeff<<std::endl;\
00319     }\
00320     MPI_Barrier(comm);\
00321   } /*end if active*/\
00322   da->vecRestoreBuffer<double>((*matPropVecPtr), matPropArr,\
00323       true, true, false, 2);\
00324 }
00325 
00326 #define FINE_TO_COARSE_BLOCK {\
00327   damg_i = damg[i];\
00328   damg_i->user = ctx;\
00329   da = damg_i->da;\
00330   assert(da->iAmActive() == daf->iAmActive());\
00331   da->createVector<double>((*matPropVecPtr), true, true, 2);\
00332   for(unsigned int j = 0; j < matPropVecPtr->size(); j++) {\
00333     (*matPropVecPtr)[j] = 0.0;\
00334   }\
00335   comm = da->getCommActive();\
00336   /*Elemental,Ghosted,Write,2 Dof.*/\
00337   matPropArr = NULL;\
00338   da->vecGetBuffer<double>((*matPropVecPtr), matPropArr,\
00339       true, true, false, 2);\
00340   double *fMatArr = NULL;\
00341   if(changedPartition) {\
00342     /*Elemental, non-Ghosted, Read-only, 2 Dof.*/\
00343     daf->vecGetBuffer<double>((*fMatPropVec), fMatArr,\
00344         true, false, true, 2);\
00345   }else {\
00346     /*Elemental, Ghosted, Read-only, 2 Dof.*/\
00347     daf->vecGetBuffer<double>((*fMatPropVec), fMatArr,\
00348         true, true, true, 2);\
00349   }\
00350   if(da->iAmActive()) {\
00351     MPI_Comm_rank(comm,&rank);\
00352     double maxCoeff = 0.0;\
00353     double minCoeff = 1.0e+8;\
00354     double globalMaxCoeff;\
00355     double globalMinCoeff;\
00356     /*Loop through the coarse and fine simultaneously*/\
00357     /*Note: If Coarse is Independent, then the*/\
00358     /*corresponding Fine is also independent.*/\
00359     /*Hence, overlapping comm with comp is possible.*/\
00360     /*First, we loop though the dependent elements.*/\
00361     /*Then we begin the communication and simulatenously*/\
00362     /*loop over the independent elements.*/\
00363     for(da->init<ot::DA_FLAGS::W_DEPENDENT>(),\
00364         daf->init<ot::DA_FLAGS::WRITABLE>();\
00365         da->curr() < da->end<ot::DA_FLAGS::W_DEPENDENT>();\
00366         da->next<ot::DA_FLAGS::W_DEPENDENT>()) {\
00367       ASSIGN_MAT_PROP_FINE_TO_COARSE_ELEM_BLOCK \
00368     } /*end dependent loop*/\
00369     da->ReadFromGhostElemsBegin<double>(matPropArr,2);\
00370     for(da->init<ot::DA_FLAGS::INDEPENDENT>(),\
00371         daf->init<ot::DA_FLAGS::WRITABLE>();\
00372         da->curr() < da->end<ot::DA_FLAGS::INDEPENDENT>();\
00373         da->next<ot::DA_FLAGS::INDEPENDENT>()) {\
00374       ASSIGN_MAT_PROP_FINE_TO_COARSE_ELEM_BLOCK \
00375     } /*end Independent loop */\
00376     da->ReadFromGhostElemsEnd<double>(matPropArr);\
00377     par::Mpi_Reduce<double>(&maxCoeff, &globalMaxCoeff, 1, MPI_MAX, 0, comm);\
00378     par::Mpi_Reduce<double>(&minCoeff, &globalMinCoeff, 1, MPI_MIN, 0, comm);\
00379     if(!rank) {\
00380       std::cout<<"Level: "<<i<<" Max Lap. Coeff: "\
00381       <<globalMaxCoeff<<" Min Lap. Coeff: "\
00382       <<globalMinCoeff<<std::endl;\
00383     }\
00384     MPI_Barrier(comm);\
00385   } /*end check if active*/\
00386   da->vecRestoreBuffer<double>((*matPropVecPtr),\
00387       matPropArr, true, true, false, 2);\
00388   if(changedPartition) {\
00389     /*Elemental, non-Ghosted, Read-only, 2 Dof.*/\
00390     daf->vecRestoreBuffer<double>((*fMatPropVec),\
00391         fMatArr, true, false, true, 2);\
00392   }else {\
00393     /*Elemental, Ghosted, Read-only, 2 Dof.*/\
00394     daf->vecRestoreBuffer<double>((*fMatPropVec),\
00395         fMatArr, true, true, true, 2);\
00396   }\
00397 }
00398 
00399 void SetUserContexts(ot::DAMG* damg) {
00400 
00401   PetscInt       jacType = 1;
00402   PetscOptionsGetInt(0, "-jacType", &jacType, 0);
00403 
00404   if(jacType == 1) { return; }
00405 
00406   PetscTruth setMatPropsAtCoarsest;
00407   PetscOptionsHasName(0,"-setMatPropsAtCoarsest",&setMatPropsAtCoarsest);
00408 
00409   if(setMatPropsAtCoarsest) {
00410     assert(jacType == 2);
00411     SetUserContextsCoarsestToFinest(damg);
00412     return;
00413   }
00414 
00415   int       nlevels = damg[0]->nlevels; //number of multigrid levels
00416 
00417   //Set Mat Props Finest to Coarsest...
00418   void * ctx;
00419 
00420   // Set for the finest level first
00421   if(jacType == 2) {
00422     ctx = new Jac2MFreeData;
00423   }else {
00424     ctx = new Jac3MFreeData;
00425   }
00426 
00427   ot::DAMG damg_i = damg[nlevels-1];
00428   int rank;
00429   MPI_Comm comm;
00430 
00431   std::vector<double> * matPropVecPtr = new std::vector<double>;
00432 
00433   if(jacType == 2) {
00434     (static_cast<Jac2MFreeData*>(ctx))->matProp = matPropVecPtr;
00435     (static_cast<Jac2MFreeData*>(ctx))->isFinestLevel = true;
00436     (static_cast<Jac2MFreeData*>(ctx))->Jmat_private = NULL;
00437     (static_cast<Jac2MFreeData*>(ctx))->inTmp = NULL;
00438     (static_cast<Jac2MFreeData*>(ctx))->outTmp = NULL;
00439   }else {
00440     (static_cast<Jac3MFreeData*>(ctx))->matProp = matPropVecPtr;
00441     (static_cast<Jac3MFreeData*>(ctx))->isFinestLevel = true;
00442     (static_cast<Jac3MFreeData*>(ctx))->isCoarsestLevel = false;
00443     (static_cast<Jac3MFreeData*>(ctx))->daf = NULL;
00444     (static_cast<Jac3MFreeData*>(ctx))->matPropFine = NULL;
00445     (static_cast<Jac3MFreeData*>(ctx))->changedPartition = false;
00446     (static_cast<Jac3MFreeData*>(ctx))->JmatThisLevel = NULL;
00447     (static_cast<Jac3MFreeData*>(ctx))->BmatThisLevel = NULL;
00448     (static_cast<Jac3MFreeData*>(ctx))->Jmat_private = NULL;
00449     (static_cast<Jac3MFreeData*>(ctx))->inTmp = NULL;
00450     (static_cast<Jac3MFreeData*>(ctx))->outTmp = NULL;
00451   }
00452 
00453   //Default values is the same as the const. coeff. case 
00454   PetscReal lapFac = 0.0;
00455   PetscReal massFac = 0.0;
00456   PetscReal massBase = 1.0;
00457   PetscTruth optFound;
00458   PetscOptionsGetReal("lap","-MatPropFac",&lapFac,&optFound);
00459   PetscOptionsGetReal("mass","-MatPropFac",&massFac,&optFound);
00460   PetscOptionsGetReal("mass","-BaseMatProp",&massBase,&optFound);
00461   unsigned int maxD ;
00462   double *matPropArr = NULL;
00463   ot::DA* da;
00464 
00465   PetscTruth setMatPropFromAnalyticFn;
00466   PetscOptionsHasName(0,"-setMatPropFromAnalyticFn",&setMatPropFromAnalyticFn);
00467 
00468   if(setMatPropFromAnalyticFn) {
00469     PetscReal lapFreq = 1.0;
00470     PetscOptionsGetReal(0,"-lapFreq",&lapFreq,&optFound);
00471     int dummyInit;
00472     SET_SINGLE_LEVEL_MAT_PROP_BLOCK(\
00473         ASSIGN_MAT_PROP_ANALYTIC_FN_ELEM_BLOCK, dummyInit = 0;) 
00474   } else {
00475     PetscTruth setCheckerBoardMatProp;
00476     PetscOptionsHasName(0,"-setCheckerBoardMatProp",&setCheckerBoardMatProp);
00477     if(setCheckerBoardMatProp) {
00478       int dummyInit;
00479       SET_SINGLE_LEVEL_MAT_PROP_BLOCK(\
00480           ASSIGN_MAT_PROP_CHECKER_BOARD_ELEM_BLOCK, dummyInit = 0;) 
00481     } else {
00482       PetscInt numCubes = 1;
00483       PetscOptionsGetInt(0,"-numCubes",&numCubes,0);
00484       if(numCubes == 1) {
00485         int dummyInit;
00486         SET_SINGLE_LEVEL_MAT_PROP_BLOCK(\
00487             ASSIGN_MAT_PROP_1_CUBE_ELEM_BLOCK, dummyInit = 0;) 
00488       } else {
00489         int dummyInit;
00490         SET_SINGLE_LEVEL_MAT_PROP_BLOCK(\
00491             ASSIGN_MAT_PROP_MULTIPLE_CUBES_ELEM_BLOCK, dummyInit = 0;) 
00492       }
00493     }
00494   }
00495 
00496   ot::DA* daf;
00497   std::vector<double>* fMatPropVec = NULL;
00498   bool changedPartition;
00499 
00500   if(nlevels > 1) {
00501     CHK_AND_SCATTER_FINE_TO_COARSE_BLOCK 
00502   }
00503 
00504   //Coarser levels
00505   for(int i = (nlevels-2); i >= 0; i--) {
00506     if(jacType == 2) {
00507       ctx = new Jac2MFreeData;
00508     }else {
00509       ctx = new Jac3MFreeData;
00510     }
00511 
00512     matPropVecPtr = new std::vector<double>;
00513 
00514     if(jacType == 2) {
00515       (static_cast<Jac2MFreeData*>(ctx))->matProp = matPropVecPtr;
00516       (static_cast<Jac2MFreeData*>(ctx))->isFinestLevel = false;
00517       (static_cast<Jac2MFreeData*>(ctx))->Jmat_private = NULL;
00518       (static_cast<Jac2MFreeData*>(ctx))->inTmp = NULL;
00519       (static_cast<Jac2MFreeData*>(ctx))->outTmp = NULL;
00520     }else {
00521       (static_cast<Jac3MFreeData*>(ctx))->matProp = matPropVecPtr;
00522       (static_cast<Jac3MFreeData*>(ctx))->isFinestLevel = false;
00523       (static_cast<Jac3MFreeData*>(ctx))->isCoarsestLevel = false;
00524       (static_cast<Jac3MFreeData*>(ctx))->daf = daf;
00525       //Note, while using matPropFine, first check changedPartition. If the
00526       //partition was changed, then the vector is a non-ghosted elemental
00527       //vector. Else it is a ghosted elemental vector.
00528       (static_cast<Jac3MFreeData*>(ctx))->matPropFine = fMatPropVec;
00529       //If the partition was changed a new vector would have been created. Else
00530       // only the pointer is copied.
00531       (static_cast<Jac3MFreeData*>(ctx))->changedPartition = changedPartition;
00532       (static_cast<Jac3MFreeData*>(ctx))->JmatThisLevel = NULL;
00533       (static_cast<Jac3MFreeData*>(ctx))->BmatThisLevel = NULL;
00534       (static_cast<Jac3MFreeData*>(ctx))->Jmat_private = NULL;
00535       (static_cast<Jac3MFreeData*>(ctx))->inTmp = NULL;
00536       (static_cast<Jac3MFreeData*>(ctx))->outTmp = NULL;
00537     }
00538 
00539     FINE_TO_COARSE_BLOCK
00540 
00541       if(changedPartition && (jacType == 2)) {
00542         fMatPropVec->clear();
00543         delete fMatPropVec;
00544       }
00545 
00546     if(i) {     
00547       CHK_AND_SCATTER_FINE_TO_COARSE_BLOCK 
00548     }
00549   }//end for i
00550 
00551   //To handle the case of a single grid, it is best to set the coarsest level's
00552   //flag separately at the end.
00553   if(jacType == 3) {
00554     (static_cast<Jac3MFreeData*>(damg[0]->user))->isCoarsestLevel = true;
00555   }
00556 }//end fn.
00557 
00558 void SetUserContextsFromPts(ot::DAMG* damg,
00559     const std::vector<double>& pts,
00560     const std::vector<double> & lapJump) {
00561 
00562   PetscInt       jacType = 1;
00563   PetscOptionsGetInt(0, "-jacType", &jacType, 0);
00564 
00565   assert(jacType != 1);
00566 
00567   assert(pts.size() == (3*lapJump.size()));
00568 
00569   int nlevels = damg[0]->nlevels; //number of mg levels
00570   ot::DAMG damg_i = damg[nlevels-1];
00571 
00572   //Set Mat Props Finest to Coarsest...
00573   void * ctx;
00574 
00575   // Set for the finest level first
00576   if(jacType == 2) {
00577     ctx = new Jac2MFreeData;
00578   }else {
00579     ctx = new Jac3MFreeData;
00580   }
00581 
00582   std::vector<double> * matPropVecPtr = new std::vector<double>;
00583 
00584   if(jacType == 2) {
00585     (static_cast<Jac2MFreeData*>(ctx))->matProp = matPropVecPtr;
00586     (static_cast<Jac2MFreeData*>(ctx))->isFinestLevel = true;
00587     (static_cast<Jac2MFreeData*>(ctx))->Jmat_private = NULL;
00588     (static_cast<Jac2MFreeData*>(ctx))->inTmp = NULL;
00589     (static_cast<Jac2MFreeData*>(ctx))->outTmp = NULL;
00590   }else {
00591     (static_cast<Jac3MFreeData*>(ctx))->matProp = matPropVecPtr;
00592     (static_cast<Jac3MFreeData*>(ctx))->isFinestLevel = true;
00593     (static_cast<Jac3MFreeData*>(ctx))->isCoarsestLevel = false;
00594     (static_cast<Jac3MFreeData*>(ctx))->daf = NULL;
00595     (static_cast<Jac3MFreeData*>(ctx))->matPropFine = NULL;
00596     (static_cast<Jac3MFreeData*>(ctx))->changedPartition = false;
00597     (static_cast<Jac3MFreeData*>(ctx))->JmatThisLevel = NULL;
00598     (static_cast<Jac3MFreeData*>(ctx))->BmatThisLevel = NULL;
00599     (static_cast<Jac3MFreeData*>(ctx))->Jmat_private = NULL;
00600     (static_cast<Jac3MFreeData*>(ctx))->inTmp = NULL;
00601     (static_cast<Jac3MFreeData*>(ctx))->outTmp = NULL;
00602   }
00603 
00604   //Default values is the same as the const. coeff. case 
00605   PetscReal lapBase = 1.0;
00606   PetscReal massBase = 1.0;
00607   PetscTruth optFound;
00608   PetscOptionsGetReal("lap","-BaseMatProp",&lapBase,&optFound);
00609   PetscOptionsGetReal("mass","-BaseMatProp",&massBase,&optFound);
00610   unsigned int maxD ;
00611   double *matPropArr = NULL;
00612   ot::DA* da;
00613   int rank;
00614   MPI_Comm comm;
00615 
00616   unsigned int ptsCtr;
00617   SET_SINGLE_LEVEL_MAT_PROP_BLOCK(ASSIGN_MAT_PROP_FROM_PTS_ELEM_BLOCK,
00618       ptsCtr = 0;) 
00619 
00620     //The coarser levels...
00621     ot::DA* daf;
00622   std::vector<double>* fMatPropVec = NULL;
00623   bool changedPartition;
00624 
00625   if(nlevels > 1) {
00626     CHK_AND_SCATTER_FINE_TO_COARSE_BLOCK 
00627   }
00628 
00629   //Coarser levels
00630   for(int i = (nlevels-2); i >= 0; i--) {
00631     if(jacType == 2) {
00632       ctx = new Jac2MFreeData;
00633     }else {
00634       ctx = new Jac3MFreeData;
00635     }
00636 
00637     matPropVecPtr = new std::vector<double>;
00638 
00639     if(jacType == 2) {
00640       (static_cast<Jac2MFreeData*>(ctx))->matProp = matPropVecPtr;
00641       (static_cast<Jac2MFreeData*>(ctx))->isFinestLevel = false;
00642       (static_cast<Jac2MFreeData*>(ctx))->Jmat_private = NULL;
00643       (static_cast<Jac2MFreeData*>(ctx))->inTmp = NULL;
00644       (static_cast<Jac2MFreeData*>(ctx))->outTmp = NULL;
00645     }else {
00646       (static_cast<Jac3MFreeData*>(ctx))->matProp = matPropVecPtr;
00647       (static_cast<Jac3MFreeData*>(ctx))->isFinestLevel = false;
00648       (static_cast<Jac3MFreeData*>(ctx))->isCoarsestLevel = false;
00649       (static_cast<Jac3MFreeData*>(ctx))->daf = daf;
00650       //Note, while using matPropFine, first check changedPartition. If the
00651       //partition was changed, then the vector is a non-ghosted elemental
00652       //vector. Else it is a ghosted elemental vector.
00653       (static_cast<Jac3MFreeData*>(ctx))->matPropFine = fMatPropVec;
00654       //If the partition was changed a new vector would have been created. Else
00655       // only the pointer is copied.
00656       (static_cast<Jac3MFreeData*>(ctx))->changedPartition = changedPartition;
00657       (static_cast<Jac3MFreeData*>(ctx))->JmatThisLevel = NULL;
00658       (static_cast<Jac3MFreeData*>(ctx))->BmatThisLevel = NULL;
00659       (static_cast<Jac3MFreeData*>(ctx))->Jmat_private = NULL;
00660       (static_cast<Jac3MFreeData*>(ctx))->inTmp = NULL;
00661       (static_cast<Jac3MFreeData*>(ctx))->outTmp = NULL;
00662     }
00663 
00664     FINE_TO_COARSE_BLOCK
00665 
00666       if(changedPartition && (jacType == 2)) {
00667         fMatPropVec->clear();
00668         delete fMatPropVec;
00669       }
00670 
00671     if(i) {     
00672       CHK_AND_SCATTER_FINE_TO_COARSE_BLOCK 
00673     }
00674   }//end for i
00675 
00676   //To handle the case of a single grid, it is best to set the coarsest level's
00677   //flag separately at the end.
00678   if(jacType == 3) {
00679     (static_cast<Jac3MFreeData*>(damg[0]->user))->isCoarsestLevel = true;
00680   }
00681 }//end fn.
00682 
00683 #undef ASSIGN_MAT_PROP_FINE_TO_COARSE_ELEM_BLOCK 
00684 #undef CHK_AND_SCATTER_FINE_TO_COARSE_BLOCK 
00685 #undef FINE_TO_COARSE_BLOCK
00686 
00687 #define ASSIGN_MAT_PROP_COARSE_TO_FINE_ELEM_BLOCK {\
00688   /*The fine loop is always Writable, but the coarse loop*/\
00689   /*could be Independent or W_Dependent. Hence the fine counter must*/\
00690   /*be incremented properly to align with the coarse.*/\
00691   unsigned int idxC= dac->curr();\
00692   Point Cpt = dac->getCurrentOffset();\
00693   assert(daf->curr() < daf->end<ot::DA_FLAGS::WRITABLE>());\
00694   while(daf->getCurrentOffset() != Cpt) {\
00695     daf->next<ot::DA_FLAGS::WRITABLE>();\
00696     assert(daf->curr() < daf->end<ot::DA_FLAGS::WRITABLE>());\
00697   }\
00698   if(daf->getLevel(daf->curr()) == dac->getLevel(idxC)) {\
00699     /*The coarse and fine elements are the same,*/\
00700     assert(daf->curr() < daf->end<ot::DA_FLAGS::WRITABLE>());\
00701     unsigned int idxF = daf->curr();\
00702     fMatArr[2*idxF] = cMatArr[2*idxC];\
00703     fMatArr[(2*idxF)+1] = cMatArr[(2*idxC)+1];\
00704     if(fMatArr[2*idxF] > maxCoeff) {\
00705       maxCoeff = fMatArr[2*idxF];\
00706     }\
00707     if(fMatArr[2*idxF] < minCoeff) {\
00708       minCoeff = fMatArr[2*idxF];\
00709     }\
00710     daf->next<ot::DA_FLAGS::WRITABLE>();\
00711   }else {\
00712     for(unsigned char cNumFine = 0; cNumFine < 8; cNumFine++) {\
00713       /*The coarse and fine elements are NOT the same. */\
00714       /*Loop over each of the 8 children of the coarse element.*/\
00715       /*These are the underlying fine elements.*/\
00716       assert(daf->curr() < daf->end<ot::DA_FLAGS::WRITABLE>());\
00717       unsigned int idxF = daf->curr();\
00718       fMatArr[2*idxF] = cMatArr[2*idxC];\
00719       fMatArr[(2*idxF)+1] = cMatArr[(2*idxC)+1];\
00720       if(fMatArr[2*idxF] > maxCoeff) {\
00721         maxCoeff = fMatArr[2*idxF];\
00722       }\
00723       if(fMatArr[2*idxF] < minCoeff) {\
00724         minCoeff = fMatArr[2*idxF];\
00725       }\
00726       daf->next<ot::DA_FLAGS::WRITABLE>();\
00727     }\
00728   }\
00729 }
00730 
00731 #define COARSE_TO_FINE_BLOCK {\
00732   /*The finer levels... */\
00733   for(int i = 1; i < nlevels; i++) {\
00734     ot::DA* dac = damg[i-1]->da;\
00735     ot::DA* daf;\
00736     bool changedPartition;\
00737     if(damg[i]->da_aux) {\
00738       daf = damg[i]->da_aux;\
00739       changedPartition = true;\
00740     }else {\
00741       daf = damg[i]->da;\
00742       changedPartition = false;\
00743     }\
00744     assert(dac->iAmActive() == daf->iAmActive());\
00745     comm = daf->getCommActive();\
00746     if(daf->iAmActive()) {\
00747       MPI_Comm_rank(comm,&rank);\
00748     }\
00749     std::vector<double>* cMatPropVec = \
00750     (static_cast<Jac2MFreeData*>(damg[i-1]->user))->matProp;\
00751     std::vector<double>* fMatPropVec = new std::vector<double>;\
00752     if(changedPartition) {\
00753       /*Elemental and Non-Ghosted*/\
00754       daf->createVector<double>((*fMatPropVec), true, false, 2);\
00755     } else {\
00756       /*Elemental and Ghosted*/\
00757       daf->createVector<double>((*fMatPropVec), true, true, 2);\
00758     }\
00759     for(unsigned int j = 0; j < fMatPropVec->size(); j++) {\
00760       (*fMatPropVec)[j] = 0.0;\
00761     }\
00762     double *cMatArr = NULL;\
00763     double *fMatArr = NULL;\
00764     double maxCoeff = 0.0;\
00765     double minCoeff = 1.0e+8;\
00766     double globalMaxCoeff;\
00767     double globalMinCoeff;\
00768     /*Read-only buffer*/\
00769     dac->vecGetBuffer<double>((*cMatPropVec), cMatArr,true, true, true, 2);\
00770     if(changedPartition) {\
00771       if(daf->iAmActive()) {\
00772         /*Writable buffer*/\
00773         daf->vecGetBuffer<double>((*fMatPropVec), fMatArr,\
00774             true, false, false, 2);\
00775         /*Can not overlap comm and comp here. So direct WRITABLE loop*/\
00776         for(dac->init<ot::DA_FLAGS::WRITABLE>(),\
00777             daf->init<ot::DA_FLAGS::WRITABLE>();\
00778             dac->curr() < dac->end<ot::DA_FLAGS::WRITABLE>();\
00779             dac->next<ot::DA_FLAGS::WRITABLE>()) {\
00780           ASSIGN_MAT_PROP_COARSE_TO_FINE_ELEM_BLOCK \
00781         } /*end loop*/\
00782         daf->vecRestoreBuffer<double>((*fMatPropVec), fMatArr,\
00783             true, false, false, 2);\
00784       }/*end if active*/\
00785       std::vector<double> tmpVecForScatter;\
00786       /*Scatter from fMatPropVec to*/\
00787       /*tmpVecForScatter (created in the function)*/\
00788       par::scatterValues<double>((*fMatPropVec), tmpVecForScatter,\
00789           (2*(damg[i]->da->getElementSize())),\
00790           daf->getComm());\
00791       fMatPropVec->clear();\
00792       delete fMatPropVec;\
00793       std::vector<double> *tmpFineVec = new std::vector<double>;\
00794       /*Elemental and Ghosted*/\
00795       damg[i]->da->createVector<double>((*tmpFineVec), true, true, 2);\
00796       if(damg[i]->da->iAmActive()) {\
00797         /*Writable buffer*/\
00798         damg[i]->da->vecGetBuffer<double>((*tmpFineVec), fMatArr,\
00799             true, true, false, 2);\
00800         /*Read-only buffer*/\
00801         double *tmpFmatArr = NULL;\
00802         damg[i]->da->vecGetBuffer<double>(tmpVecForScatter, tmpFmatArr,\
00803             true, false, true, 2);\
00804         /*Copy from tmpVecForScatter to tmpFineVec*/\
00805         /*W_DEPENDENT loop*/\
00806         for(damg[i]->da->init<ot::DA_FLAGS::W_DEPENDENT>();\
00807             damg[i]->da->curr() < damg[i]->da->end<ot::DA_FLAGS::W_DEPENDENT>();\
00808             damg[i]->da->next<ot::DA_FLAGS::W_DEPENDENT>()) {\
00809           unsigned int idxCurr = damg[i]->da->curr();\
00810           fMatArr[2*idxCurr] = tmpFmatArr[2*idxCurr]; \
00811           fMatArr[(2*idxCurr) + 1] = tmpFmatArr[(2*idxCurr) + 1]; \
00812         } /*end dependent loop*/\
00813         /*Begin read from ghost elements on the da grid*/      \
00814         damg[i]->da->ReadFromGhostElemsBegin<double>(fMatArr,2);\
00815         /*Overlap communication with Independent loop*/\
00816         for(damg[i]->da->init<ot::DA_FLAGS::INDEPENDENT>();\
00817             damg[i]->da->curr() < damg[i]->da->end<ot::DA_FLAGS::INDEPENDENT>();\
00818             damg[i]->da->next<ot::DA_FLAGS::INDEPENDENT>()) {\
00819           unsigned int idxCurr = damg[i]->da->curr();\
00820           fMatArr[2*idxCurr] = tmpFmatArr[2*idxCurr]; \
00821           fMatArr[(2*idxCurr) + 1] = tmpFmatArr[(2*idxCurr) + 1]; \
00822         } /*end Independent loop */\
00823         /*End read from ghost elements on the da grid*/      \
00824         damg[i]->da->ReadFromGhostElemsEnd<double>(fMatArr);\
00825         damg[i]->da->vecRestoreBuffer<double>(tmpVecForScatter, tmpFmatArr,\
00826             true, false, true, 2);\
00827         tmpVecForScatter.clear();\
00828         damg[i]->da->vecRestoreBuffer<double>((*tmpFineVec), fMatArr,\
00829             true, true, false, 2);\
00830       }/*end if active*/\
00831       fMatPropVec = tmpFineVec;\
00832       tmpFineVec = NULL;      \
00833     } else {\
00834       if(daf->iAmActive()) {\
00835         /*Writable buffer*/\
00836         daf->vecGetBuffer<double>((*fMatPropVec), fMatArr,\
00837             true, true, false, 2);\
00838         /*W_DEPENDENT loop*/\
00839         for(dac->init<ot::DA_FLAGS::W_DEPENDENT>(),\
00840             daf->init<ot::DA_FLAGS::WRITABLE>();\
00841             dac->curr() < dac->end<ot::DA_FLAGS::W_DEPENDENT>();\
00842             dac->next<ot::DA_FLAGS::W_DEPENDENT>()) {\
00843           ASSIGN_MAT_PROP_COARSE_TO_FINE_ELEM_BLOCK \
00844         } /*end dependent loop*/\
00845         daf->ReadFromGhostElemsBegin<double>(fMatArr,2);\
00846         /*Overlap communication with INDEPEDENT*/\
00847         for(dac->init<ot::DA_FLAGS::INDEPENDENT>(),\
00848             daf->init<ot::DA_FLAGS::WRITABLE>();\
00849             dac->curr() < dac->end<ot::DA_FLAGS::INDEPENDENT>();\
00850             dac->next<ot::DA_FLAGS::INDEPENDENT>()) {\
00851           ASSIGN_MAT_PROP_COARSE_TO_FINE_ELEM_BLOCK \
00852         } /*end Independent loop */\
00853         daf->ReadFromGhostElemsEnd<double>(fMatArr);\
00854         daf->vecRestoreBuffer<double>((*fMatPropVec), fMatArr,\
00855             true, true, false, 2);\
00856       } /*end check if active*/\
00857     }/*end if changedPartition*/\
00858     dac->vecRestoreBuffer<double>((*cMatPropVec), cMatArr,\
00859         true, true, true, 2);\
00860     ctx = new Jac2MFreeData;\
00861     (static_cast<Jac2MFreeData*>(ctx))->matProp = fMatPropVec;\
00862     (static_cast<Jac2MFreeData*>(ctx))->isFinestLevel = (i == (nlevels-1));\
00863     (static_cast<Jac2MFreeData*>(ctx))->Jmat_private = NULL;\
00864     (static_cast<Jac2MFreeData*>(ctx))->inTmp = NULL;\
00865     (static_cast<Jac2MFreeData*>(ctx))->outTmp = NULL;\
00866     damg[i]->user = ctx;\
00867     if(daf->iAmActive()) {\
00868       par::Mpi_Reduce<double>(&maxCoeff, &globalMaxCoeff, 1, MPI_MAX, 0, comm);\
00869       par::Mpi_Reduce<double>(&minCoeff, &globalMinCoeff, 1, MPI_MIN, 0, comm);\
00870       if(!rank) {\
00871         std::cout<<"Level: "<<i<<" Max Lap. Coeff: "\
00872         <<globalMaxCoeff<<" Min Lap. Coeff: "\
00873         <<globalMinCoeff<<std::endl;\
00874       }\
00875       MPI_Barrier(comm);\
00876     }/*end if active*/\
00877   }/*end for i*/\
00878 }
00879 
00880 void SetUserContextsCoarsestToFinest(ot::DAMG* damg) {
00881   //jactype =2 only
00882   int       nlevels = damg[0]->nlevels; //number of multigrid levels
00883 
00884   //coarsest first...
00885   void * ctx = new Jac2MFreeData;   
00886   ot::DAMG damg_i = damg[0];
00887 
00888   std::vector<double> * matPropVecPtr = new std::vector<double>;
00889   (static_cast<Jac2MFreeData*>(ctx))->matProp = matPropVecPtr;
00890   (static_cast<Jac2MFreeData*>(ctx))->isFinestLevel = (nlevels == 1);
00891   (static_cast<Jac2MFreeData*>(ctx))->Jmat_private = NULL;
00892   (static_cast<Jac2MFreeData*>(ctx))->inTmp = NULL;
00893   (static_cast<Jac2MFreeData*>(ctx))->outTmp = NULL;
00894 
00895   int rank;
00896   MPI_Comm comm;
00897 
00898   //Default values is the same as the const. coeff. case 
00899   PetscReal lapFac = 0.0;
00900   PetscReal massFac = 0.0;
00901   PetscReal massBase = 1.0;
00902   PetscTruth optFound;
00903   PetscOptionsGetReal("lap","-MatPropFac",&lapFac,&optFound);
00904   PetscOptionsGetReal("mass","-MatPropFac",&massFac,&optFound);
00905   PetscOptionsGetReal("mass","-BaseMatProp",&massBase,&optFound);
00906   unsigned int maxD;
00907   ot::DA* da;
00908   double *matPropArr = NULL;
00909 
00910   PetscTruth setMatPropFromAnalyticFn;
00911   PetscOptionsHasName(0,"-setMatPropFromAnalyticFn",&setMatPropFromAnalyticFn);
00912 
00913   if(setMatPropFromAnalyticFn) {
00914     PetscReal lapFreq = 1.0;
00915     PetscOptionsGetReal(0,"-lapFreq",&lapFreq,&optFound);
00916     int dummyInit;
00917     SET_SINGLE_LEVEL_MAT_PROP_BLOCK(\
00918         ASSIGN_MAT_PROP_ANALYTIC_FN_ELEM_BLOCK, dummyInit = 0;) 
00919   } else {
00920     PetscTruth setCheckerBoardMatProp;
00921     PetscOptionsHasName(0,"-setCheckerBoardMatProp",&setCheckerBoardMatProp);
00922     if(setCheckerBoardMatProp) {
00923       int dummyInit;
00924       SET_SINGLE_LEVEL_MAT_PROP_BLOCK(ASSIGN_MAT_PROP_CHECKER_BOARD_ELEM_BLOCK, 
00925           dummyInit = 0;) 
00926     } else {
00927       PetscInt numCubes = 1;
00928       PetscOptionsGetInt(0,"-numCubes",&numCubes,0);
00929       if(numCubes == 1) {
00930         int dummyInit;
00931         SET_SINGLE_LEVEL_MAT_PROP_BLOCK(ASSIGN_MAT_PROP_1_CUBE_ELEM_BLOCK, 
00932             dummyInit = 0;) 
00933       } else {
00934         int dummyInit;
00935         SET_SINGLE_LEVEL_MAT_PROP_BLOCK(ASSIGN_MAT_PROP_MULTIPLE_CUBES_ELEM_BLOCK,
00936             dummyInit = 0;) 
00937       }
00938     }
00939   }
00940 
00941   COARSE_TO_FINE_BLOCK
00942 
00943 }//end fn.
00944 
00945 void SetCoarseToFineFromPts(ot::DAMG* damg,
00946     const std::vector<double>& pts,
00947     const std::vector<double> & lapJump) {
00948 
00949   assert(pts.size() == (3*lapJump.size()));
00950 
00951   //jactype =2 only
00952   int       nlevels = damg[0]->nlevels; //number of multigrid levels
00953 
00954   //coarsest first...
00955   void * ctx = new Jac2MFreeData;   
00956   ot::DAMG damg_i = damg[0];
00957 
00958   std::vector<double> * matPropVecPtr = new std::vector<double>;
00959   (static_cast<Jac2MFreeData*>(ctx))->matProp = matPropVecPtr;
00960   (static_cast<Jac2MFreeData*>(ctx))->isFinestLevel = (nlevels == 1);
00961   (static_cast<Jac2MFreeData*>(ctx))->Jmat_private = NULL;
00962   (static_cast<Jac2MFreeData*>(ctx))->inTmp = NULL;
00963   (static_cast<Jac2MFreeData*>(ctx))->outTmp = NULL;
00964 
00965   int rank;
00966   MPI_Comm comm;
00967 
00968   //Default values is the same as the const. coeff. case 
00969   PetscReal lapBase = 1.0;
00970   PetscReal massBase = 1.0;
00971   PetscTruth optFound;
00972   PetscOptionsGetReal("lap","-BaseMatProp",&lapBase,&optFound);
00973   PetscOptionsGetReal("mass","-BaseMatProp",&massBase,&optFound);
00974   unsigned int maxD;
00975   ot::DA* da;
00976   double *matPropArr = NULL;
00977 
00978   unsigned int ptsCtr;
00979   SET_SINGLE_LEVEL_MAT_PROP_BLOCK(ASSIGN_MAT_PROP_FROM_PTS_ELEM_BLOCK,
00980       ptsCtr = 0;) 
00981 
00982     COARSE_TO_FINE_BLOCK
00983 
00984 }//end fn.
00985 
00986 #undef ASSIGN_MAT_PROP_FROM_PTS_ELEM_BLOCK 
00987 #undef ASSIGN_MAT_PROP_1_CUBE_ELEM_BLOCK 
00988 #undef ASSIGN_MAT_PROP_MULTIPLE_CUBES_ELEM_BLOCK 
00989 #undef ASSIGN_MAT_PROP_CHECKER_BOARD_ELEM_BLOCK 
00990 
00991 #undef SET_SINGLE_LEVEL_MAT_PROP_BLOCK
00992 
00993 #undef ASSIGN_MAT_PROP_COARSE_TO_FINE_ELEM_BLOCK 
00994 #undef COARSE_TO_FINE_BLOCK 
00995 
00996 #undef ASSIGN_MAT_PROP_ANALYTIC_FN_ELEM_BLOCK 
00997 #undef SQUARE
00998 
00999 void DestroyUserContexts(ot::DAMG* damg) {
01000 
01001   PetscInt       jacType = 1;
01002   PetscOptionsGetInt(0,"-jacType",&jacType,0);
01003 
01004   if(jacType == 1) { return; }
01005 
01006   int       nlevels = damg[0]->nlevels; //number of multigrid levels
01007 
01008   if(jacType == 2) {
01009     for(int i = 0; i < nlevels; i++) {
01010       Jac2MFreeData* ctx = (static_cast<Jac2MFreeData*>(damg[i]->user));
01011       ctx->matProp->clear();
01012       delete ctx->matProp;
01013       ctx->matProp = NULL;
01014       if(ctx->Jmat_private) {
01015         MatDestroy(ctx->Jmat_private);
01016         ctx->Jmat_private = NULL;
01017       }
01018       if(ctx->inTmp) {
01019         VecDestroy(ctx->inTmp);
01020         ctx->inTmp = NULL;
01021       }
01022       if(ctx->outTmp) {
01023         VecDestroy(ctx->outTmp);
01024         ctx->outTmp = NULL;
01025       }
01026       delete ctx;
01027       ctx = NULL;
01028     }
01029   }//end type-2
01030 
01031   if(jacType == 3) {
01032     for(int i = 0; i < nlevels; i++) {
01033       Jac3MFreeData* ctx = (static_cast<Jac3MFreeData*>(damg[i]->user));
01034       //The coarsest level will not create matProp. It will only use the finer
01035       //level's material properties.
01036       if(ctx->matProp) {
01037         ctx->matProp->clear();
01038         delete ctx->matProp;
01039         ctx->matProp = NULL;
01040       }
01041       //Sometimes, memory is allocated for the fine material property vector.
01042       //Sometimes it is just a copy of pointers.
01043       if(ctx->changedPartition) {
01044         ctx->matPropFine->clear();
01045         delete ctx->matPropFine;
01046         ctx->matPropFine = NULL;
01047       }
01048       if(ctx->Jmat_private) {
01049         MatDestroy(ctx->Jmat_private);
01050         ctx->Jmat_private = NULL;
01051       }
01052       if(ctx->inTmp) {
01053         VecDestroy(ctx->inTmp);
01054         ctx->inTmp = NULL;
01055       }
01056       if(ctx->outTmp) {
01057         VecDestroy(ctx->outTmp);
01058         ctx->outTmp = NULL;
01059       }
01060       delete ctx;
01061       ctx = NULL;
01062     }
01063   }//end type-3
01064 
01065 }//end fn.
01066 
01067 
01068