Mesh.cc 34.5 KB
Newer Older
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
    // === Write periodic associations. ===
    int mapSize = periodicAssociations.size();
    SerUtil::serialize(out, mapSize);
    for (std::map<BoundaryType, VertexVector*>::iterator it = periodicAssociations.begin();
	 it != periodicAssociations.end(); ++it) {
      BoundaryType b = it->first;

      // Check which DOFAdmin is used for the current VertexVector we want to serialize.
      int ithAdmin = -1;
      for (int i = 0; i < static_cast<int>(admin.size()); i++) {
	if (admin[i] == it->second->getAdmin()) {
	  ithAdmin = i;
	  break;
	}
      }
      TEST_EXIT(ithAdmin >= 0)
	("No DOFAdmin found for serialization of periodic associations!\n");

      SerUtil::serialize(out, b);
      SerUtil::serialize(out, ithAdmin);
      it->second->serialize(out);
    }

1024
1025
1026
    serializedDOFs.clear();
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1027
  void Mesh::deserialize(std::istream &in)
1028
  {
1029
1030
    FUNCNAME("Mesh::deserialize()");

1031
1032
1033
1034
1035
1036
    serializedDOFs.clear();

    in >> name;
    in.get();

    int oldVal = dim;
1037
1038
    SerUtil::deserialize(in, dim);
    TEST_EXIT_DBG(oldVal == 0 || dim == oldVal)("Invalid dimension!\n");
1039

1040
1041
1042
1043
1044
1045
    SerUtil::deserialize(in, nVertices);
    SerUtil::deserialize(in, nEdges);
    SerUtil::deserialize(in, nLeaves);
    SerUtil::deserialize(in, nElements);
    SerUtil::deserialize(in, nFaces);
    SerUtil::deserialize(in, maxEdgeNeigh);
1046
1047
1048

    diam.deserialize(in);

1049
    SerUtil::deserialize(in, preserveCoarseDOFs);
1050
1051

    oldVal = nDOFEl;
1052
1053
    SerUtil::deserialize(in, nDOFEl);
    TEST_EXIT_DBG(oldVal == 0 || nDOFEl == oldVal)("Invalid nDOFEl!\n");
1054
1055
1056
1057

    nDOF.deserialize(in);

    oldVal = nNodeEl;
1058
1059
    SerUtil::deserialize(in, nNodeEl);
    TEST_EXIT_DBG(oldVal == 0 || nNodeEl == oldVal)("Invalid nNodeEl!\n");
1060
1061
1062

    node.deserialize(in);

1063
1064
1065

    // === Read admins. ===

1066
    int size;
1067
    SerUtil::deserialize(in, size);
1068
    admin.resize(size, NULL);
1069
    for (int i = 0; i < size; i++) {
1070
      if (!admin[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1071
	admin[i] = new DOFAdmin(this);
1072

1073
1074
1075
      admin[i]->deserialize(in);
    }

1076
    SerUtil::deserialize(in, size);
Thomas Witkowski's avatar
Thomas Witkowski committed
1077
    std::vector< std::vector<int> > neighbourIndices(size);
1078

1079
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++)
1080
      if (macroElements[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1081
	delete macroElements[i];
1082
1083
1084
1085
1086
1087
1088
1089

    // All macro elements are stored in the list \ref macroElements with continous 
    // index from 0 to n - 1. But the macro element index may be different and may
    // contain holes (for example when macro elements were removed because of domain
    // decomposition based parallelization. Therefore we create a temporary map
    // from macro element indices to the continous index of \ref macroElements. This
    // will be used later to find the correct neighbours of the macro elements.
    std::map<int, int> elIndexVecIndex;
1090

1091
    macroElements.resize(size);
1092
    for (int i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1093
      macroElements[i] = new MacroElement(dim);
1094
1095
      macroElements[i]->writeNeighboursTo(&(neighbourIndices[i]));
      macroElements[i]->deserialize(in);
1096
      elIndexVecIndex[macroElements[i]->getIndex()] = i;
1097
1098
    }

1099
1100
    SerUtil::deserialize(in, elementIndex);
    SerUtil::deserialize(in, initialized);
1101
1102

    // set neighbour pointer in macro elements
1103
    int nNeighbour = getGeo(NEIGH);
1104
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++) {
1105
      for (int j = 0; j < nNeighbour; j++) {
1106
	int index = neighbourIndices[i][j];
1107
1108
1109
1110
1111
1112
1113

	if (index != -1) {
	  TEST_EXIT_DBG(elIndexVecIndex.count(index) == 1)
	    ("Cannot find correct index from neighbouring macro element!\n");

	  index = elIndexVecIndex[index];

1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
	  macroElements[i]->setNeighbour(j, macroElements[index]);
	} else {
	  macroElements[i]->setNeighbour(j, NULL);
	}
      }
    }

    // set mesh pointer in elements
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(this, -1, CALL_EVERY_EL_PREORDER);
1124
    while (elInfo) {
1125
1126
1127
1128
1129
      elInfo->getElement()->setMesh(this);
      elInfo = stack.traverseNext(elInfo);
    }

    serializedDOFs.clear();
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145

    
    /// === Read periodic assoications. ===

    int mapSize = 0;
    SerUtil::deserialize(in, mapSize);
    for (int i = 0; i < mapSize; i++) {
      BoundaryType b = 0;
      int ithAdmin = 0;
      SerUtil::deserialize(in, b);
      SerUtil::deserialize(in, ithAdmin);

      VertexVector *tmpvec = new VertexVector(admin[ithAdmin], "");
      tmpvec->deserialize(in);
      periodicAssociations[b] = tmpvec;      
    }
1146
1147
1148
1149
  }

  void Mesh::initialize() 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1150
1151
1152
    std::string macroFilename("");
    std::string valueFilename("");
    std::string periodicFile("");
1153
1154
1155
1156
1157
1158
1159
1160
1161
    int check = 1;

    GET_PARAMETER(0, name + "->macro file name",  &macroFilename);
    GET_PARAMETER(0, name + "->value file name",  &valueFilename);
    GET_PARAMETER(0, name + "->periodic file", &periodicFile);
    GET_PARAMETER(0, name + "->check", "%d", &check);
    GET_PARAMETER(0, name + "->preserve coarse dofs", "%d", &preserveCoarseDOFs);

    if (macroFilename.length()) {
1162
      macroFileInfo = MacroReader::readMacro(macroFilename.c_str(), this,
1163
1164
 					     periodicFile == "" ? NULL : periodicFile.c_str(),
 					     check);
1165
1166
1167

      // If there is no value file which should be written, we can delete
      // the information of the macro file.
1168
1169
      if (!valueFilename.length())
  	clearMacroFileInfo();
1170
1171
1172
1173
1174
    }

    initialized = true;
  }

1175
1176
  bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1177
1178
    std::map<BoundaryType, VertexVector*>::iterator it;
    std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
1179
    for (it = periodicAssociations.begin(); it != end; ++it)
1180
1181
1182
1183
1184
      if ((*(it->second))[dof1] == dof2)
	return true;
    return false;
  }

1185
1186
  bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1187
1188
1189
    std::vector<DegreeOfFreedom> associatedToDOF1;
    std::map<BoundaryType, VertexVector*>::iterator it;
    std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
1190
1191
1192
    DegreeOfFreedom dof, assDOF;

    associatedToDOF1.push_back(dof1);
Thomas Witkowski's avatar
Thomas Witkowski committed
1193
1194
1195
    for (it = periodicAssociations.begin(); it != end; ++it) {
      int size = static_cast<int>(associatedToDOF1.size());
      for (int i = 0; i < size; i++) {
1196
1197
	dof = associatedToDOF1[i];
	assDOF = (*(it->second))[dof];
Thomas Witkowski's avatar
Thomas Witkowski committed
1198
	if (assDOF == dof2) {
1199
1200
	  return true;
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1201
1202
	  if (assDOF != dof) 
	    associatedToDOF1.push_back(assDOF);
1203
1204
1205
1206
1207
1208
1209
1210
	}
      }
    }
    return false;
  }

  void Mesh::clearMacroFileInfo()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1211
    macroFileInfo->clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
1212
    delete macroFileInfo;
1213
    macroFileInfo = NULL;
1214
  }
1215

1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
  void Mesh::computeMatrixFillin(const FiniteElemSpace *feSpace, 
				 std::vector<int> &nnzInRow, int &overall, int &average)
  {
    std::map<DegreeOfFreedom, int> dofCounter;

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(this, -1, Mesh::CALL_LEAF_EL);
    ElementDofIterator elDofIter(feSpace);
    while (elInfo) {
      elDofIter.reset(elInfo->getElement());
      do {
	DegreeOfFreedom dof = elDofIter.getDof();
	if (dofCounter.count(dof) == 0) {
	  dofCounter[dof] = 1;
	} else {
	  dofCounter[dof]++;
	}
      } while (elDofIter.next());
      elInfo = stack.traverseNext(elInfo);
    } 

    overall = 0;
    for (std::map<DegreeOfFreedom, int>::iterator it = dofCounter.begin();
	 it != dofCounter.end(); ++it) {
      overall += it->second * 15;
    }
  }

1244
1245
  int Mesh::calcMemoryUsage()
  {
1246
    int result = sizeof(Mesh);
1247

1248
1249
1250
1251
1252
1253
    result += nDOFEl;
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
      result += admin[i]->calcMemoryUsage();
      result += admin[i]->getUsedSize() * sizeof(DegreeOfFreedom);
    }
    
1254
1255
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++)
      result += macroElements[i]->calcMemoryUsage();    
1256
1257
1258
    
    return result;
  }
1259
}
For faster browsing, not all history is shown. View entire blame