00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceMap.hpp"
00032 #include "PlayaTabs.hpp"
00033 #include "SundanceOut.hpp"
00034 #include "SundanceUnfoldPeriodicDF.hpp"
00035 #include "SundancePeriodicMesh1D.hpp"
00036 #include "SundancePartitionedLineMesher.hpp"
00037 #include "SundanceMeshSource.hpp"
00038 #include "SundanceBasicSimplicialMeshType.hpp"
00039 #include "SundanceDOFMapBase.hpp"
00040 #include "PlayaMPIContainerComm.hpp"
00041 #include "Teuchos_TimeMonitor.hpp"
00042 #include "PlayaLoadableVector.hpp"
00043
00044 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00045 #include "PlayaVectorSpaceImpl.hpp"
00046 #include "PlayaVectorImpl.hpp"
00047 #endif
00048
00049
00050 namespace Sundance
00051 {
00052 using namespace Teuchos;
00053 using namespace Playa;
00054
00055 Mesh unfoldPeriodicMesh(const Mesh& mesh)
00056 {
00057 const MeshBase* mb = mesh.ptr().get();
00058 const PeriodicMesh1D* pm = dynamic_cast<const PeriodicMesh1D*>(mb);
00059
00060 TEUCHOS_TEST_FOR_EXCEPT(pm==0);
00061
00062 int numElems = mesh.numCells(1);
00063 double a = mesh.nodePosition(0)[0];
00064 double b = mesh.nodePosition(numElems)[0];
00065
00066 MeshType meshType = new BasicSimplicialMeshType();
00067 int verb=0;
00068 MeshSource mesher = new PartitionedLineMesher(a, b, numElems, meshType, verb,
00069 MPIComm::self());
00070
00071 Mesh rtn = mesher.getMesh();
00072
00073 return rtn;
00074 }
00075
00076
00077 DiscreteSpace unfoldPeriodicDiscreteSpace(const DiscreteSpace& space)
00078 {
00079 VectorType<double> vecType = space.vecType();
00080 Mesh mesh = unfoldPeriodicMesh(space.mesh());
00081 BasisArray basis = space.basis();
00082
00083 return DiscreteSpace(mesh, basis, vecType);
00084 }
00085
00086 Expr unfoldPeriodicDiscreteFunction(const Expr& f, const string& name)
00087 {
00088 const DiscreteFunction* df = DiscreteFunction::discFunc(f);
00089 TEUCHOS_TEST_FOR_EXCEPT(df==0);
00090
00091
00092 DiscreteSpace perSpace = df->discreteSpace();
00093 DiscreteSpace space = unfoldPeriodicDiscreteSpace(perSpace);
00094
00095 Mesh oldMesh = perSpace.mesh();
00096 Mesh newMesh = space.mesh();
00097
00098 df->updateGhosts();
00099 RCP<GhostView<double> > ghostView = df->ghostView();
00100
00101 Vector<double> newVec = space.createVector();
00102
00103
00104 const RCP<DOFMapBase>& oldMap = perSpace.map();
00105 const RCP<DOFMapBase>& newMap = space.map();
00106
00107
00108
00109 Array<int> oldCellLID(oldMesh.numCells(1));
00110 for (int i=0; i<oldMesh.numCells(1); i++)
00111 {
00112 oldCellLID[i] = i;
00113 }
00114 RCP<const Set<int> > funcs = oldMap->allowedFuncsOnCellBatch(1, oldCellLID);
00115
00116 Array<Array<int> > oldElemDofs;
00117 Array<Array<int> > newElemDofs;
00118 Array<int> oldNNodes;
00119 RCP<const MapStructure> oldMapStruct ;
00120 RCP<const MapStructure> newMapStruct ;
00121
00122 if (funcs->size())
00123 {
00124 oldMapStruct
00125 = oldMap->getDOFsForCellBatch(1, oldCellLID, *funcs, oldElemDofs,
00126 oldNNodes, 0);
00127
00128 newMapStruct
00129 = newMap->getDOFsForCellBatch(1, oldCellLID, *funcs, newElemDofs,
00130 oldNNodes, 0);
00131
00132 for (int chunk=0; chunk<oldMapStruct->numBasisChunks(); chunk++)
00133 {
00134 int nf = oldMapStruct->numFuncs(chunk);
00135 int nDofsPerCell = oldNNodes[chunk] * nf;
00136 for (int c=0; c<oldCellLID.size(); c++)
00137 {
00138 for (int d=0; d<nDofsPerCell; d++)
00139 {
00140 int oldDof = oldElemDofs[chunk][nDofsPerCell*c + d];
00141 int newDof = newElemDofs[chunk][nDofsPerCell*c + d];
00142 loadable(newVec)->setElement(newDof, ghostView->getElement(oldDof));
00143 }
00144 }
00145 }
00146 }
00147
00148
00149 Array<int> oldVertLID(oldMesh.numCells(0));
00150 Array<int> newVertLID(newMesh.numCells(0));
00151 for (int i=0; i<oldMesh.numCells(0); i++)
00152 {
00153 oldVertLID[i] = i;
00154 }
00155 for (int i=0; i<newMesh.numCells(0); i++)
00156 {
00157 newVertLID[i] = i;
00158 }
00159 if (funcs->size())
00160 {
00161 funcs = oldMap->allowedFuncsOnCellBatch(0, oldCellLID);
00162
00163 oldMapStruct = oldMap->getDOFsForCellBatch(0, oldVertLID, *funcs,
00164 oldElemDofs,
00165 oldNNodes, 0);
00166 newMapStruct
00167 = newMap->getDOFsForCellBatch(0, newVertLID, *funcs, newElemDofs,
00168 oldNNodes, 0);
00169
00170 for (int chunk=0; chunk<oldMapStruct->numBasisChunks(); chunk++)
00171 {
00172 int nf = oldMapStruct->numFuncs(chunk);
00173 int nDofsPerCell = oldNNodes[chunk] * nf;
00174 for (int c=0; c<newVertLID.size(); c++)
00175 {
00176 if (c < oldVertLID.size())
00177 {
00178 for (int d=0; d<nDofsPerCell; d++)
00179 {
00180 int oldDof = oldElemDofs[chunk][nDofsPerCell*c + d];
00181 int newDof = newElemDofs[chunk][nDofsPerCell*c + d];
00182 loadable(newVec)->setElement(newDof, ghostView->getElement(oldDof));
00183 }
00184 }
00185 else
00186 {
00187 for (int d=0; d<nDofsPerCell; d++)
00188 {
00189 int oldDof = oldElemDofs[chunk][d];
00190 int newDof = newElemDofs[chunk][nDofsPerCell*c + d];
00191 loadable(newVec)->setElement(newDof, ghostView->getElement(oldDof));
00192 }
00193 }
00194 }
00195 }
00196 }
00197
00198 return new DiscreteFunction(space, newVec, name);
00199 }
00200
00201
00202 }
00203
00204