Mesh.cc 33.2 KB
Newer Older
1001

Thomas Witkowski's avatar
Thomas Witkowski committed
1002
  void Mesh::deserialize(std::istream &in)
1003
  {
1004
1005
    FUNCNAME("Mesh::deserialize()");

1006
1007
1008
1009
1010
1011
    serializedDOFs.clear();

    in >> name;
    in.get();

    int oldVal = dim;
1012
1013
    SerUtil::deserialize(in, dim);
    TEST_EXIT_DBG(oldVal == 0 || dim == oldVal)("Invalid dimension!\n");
1014

1015
1016
1017
1018
1019
1020
    SerUtil::deserialize(in, nVertices);
    SerUtil::deserialize(in, nEdges);
    SerUtil::deserialize(in, nLeaves);
    SerUtil::deserialize(in, nElements);
    SerUtil::deserialize(in, nFaces);
    SerUtil::deserialize(in, maxEdgeNeigh);
1021
1022
1023

    diam.deserialize(in);

1024
    SerUtil::deserialize(in, preserveCoarseDOFs);
1025
1026

    oldVal = nDOFEl;
1027
1028
    SerUtil::deserialize(in, nDOFEl);
    TEST_EXIT_DBG(oldVal == 0 || nDOFEl == oldVal)("Invalid nDOFEl!\n");
1029
1030
1031
1032

    nDOF.deserialize(in);

    oldVal = nNodeEl;
1033
1034
    SerUtil::deserialize(in, nNodeEl);
    TEST_EXIT_DBG(oldVal == 0 || nNodeEl == oldVal)("Invalid nNodeEl!\n");
1035
1036
1037

    node.deserialize(in);

1038
1039
1040

    // === Read admins. ===

1041
    int size;
1042
    SerUtil::deserialize(in, size);
1043
    admin.resize(size, NULL);
1044
    for (int i = 0; i < size; i++) {
1045
      if (!admin[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1046
	admin[i] = new DOFAdmin(this);
1047

1048
1049
1050
      admin[i]->deserialize(in);
    }

1051
    SerUtil::deserialize(in, size);
Thomas Witkowski's avatar
Thomas Witkowski committed
1052
    std::vector< std::vector<int> > neighbourIndices(size);
1053

1054
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++)
1055
      if (macroElements[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1056
	delete macroElements[i];
1057
1058
1059
1060
1061
1062
1063
1064

    // 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;
1065

1066
    macroElements.resize(size);
1067
    for (int i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1068
      macroElements[i] = new MacroElement(dim);
1069
1070
      macroElements[i]->writeNeighboursTo(&(neighbourIndices[i]));
      macroElements[i]->deserialize(in);
1071
      elIndexVecIndex[macroElements[i]->getIndex()] = i;
1072
1073
    }

1074
1075
    SerUtil::deserialize(in, elementIndex);
    SerUtil::deserialize(in, initialized);
1076
1077

    // set neighbour pointer in macro elements
1078
    int nNeighbour = getGeo(NEIGH);
1079
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++) {
1080
      for (int j = 0; j < nNeighbour; j++) {
1081
	int index = neighbourIndices[i][j];
1082
1083
1084
1085
1086
1087
1088

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

	  index = elIndexVecIndex[index];

1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
	  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);
1099
    while (elInfo) {
1100
1101
1102
1103
1104
      elInfo->getElement()->setMesh(this);
      elInfo = stack.traverseNext(elInfo);
    }

    serializedDOFs.clear();
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120

    
    /// === 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;      
    }
1121
1122
1123
1124
  }

  void Mesh::initialize() 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1125
1126
1127
    std::string macroFilename("");
    std::string valueFilename("");
    std::string periodicFile("");
1128
1129
1130
1131
1132
1133
1134
1135
1136
    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()) {
1137
1138
1139
      macroFileInfo = MacroReader::readMacro(macroFilename.c_str(), this,
					     periodicFile == "" ? NULL : periodicFile.c_str(),
					     check);
1140
1141
1142

      // If there is no value file which should be written, we can delete
      // the information of the macro file.
1143
      if (!valueFilename.length())
1144
1145
1146
1147
1148
1149
	clearMacroFileInfo();
    }

    initialized = true;
  }

1150
1151
  bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1152
1153
    std::map<BoundaryType, VertexVector*>::iterator it;
    std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
1154
    for (it = periodicAssociations.begin(); it != end; ++it)
1155
1156
1157
1158
1159
      if ((*(it->second))[dof1] == dof2)
	return true;
    return false;
  }

1160
1161
  bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1162
1163
1164
    std::vector<DegreeOfFreedom> associatedToDOF1;
    std::map<BoundaryType, VertexVector*>::iterator it;
    std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
1165
1166
1167
    DegreeOfFreedom dof, assDOF;

    associatedToDOF1.push_back(dof1);
Thomas Witkowski's avatar
Thomas Witkowski committed
1168
1169
1170
    for (it = periodicAssociations.begin(); it != end; ++it) {
      int size = static_cast<int>(associatedToDOF1.size());
      for (int i = 0; i < size; i++) {
1171
1172
	dof = associatedToDOF1[i];
	assDOF = (*(it->second))[dof];
Thomas Witkowski's avatar
Thomas Witkowski committed
1173
	if (assDOF == dof2) {
1174
1175
	  return true;
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1176
1177
	  if (assDOF != dof) 
	    associatedToDOF1.push_back(assDOF);
1178
1179
1180
1181
1182
1183
1184
1185
	}
      }
    }
    return false;
  }

  void Mesh::clearMacroFileInfo()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1186
    macroFileInfo->clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
1187
    delete macroFileInfo;
1188
    macroFileInfo = NULL;
1189
  }
1190

1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
  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;
    }
  }

1219
1220
  int Mesh::calcMemoryUsage()
  {
1221
    int result = sizeof(Mesh);
1222

1223
1224
1225
1226
1227
1228
    result += nDOFEl;
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
      result += admin[i]->calcMemoryUsage();
      result += admin[i]->getUsedSize() * sizeof(DegreeOfFreedom);
    }
    
1229
1230
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++)
      result += macroElements[i]->calcMemoryUsage();    
1231
1232
1233
    
    return result;
  }
1234
}
For faster browsing, not all history is shown. View entire blame