Mesh.cc 33.2 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1001
  void Mesh::deserialize(std::istream &in)
1002
  {
1003
1004
    FUNCNAME("Mesh::deserialize()");

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

    in >> name;
    in.get();

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

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

    diam.deserialize(in);

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

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

    nDOF.deserialize(in);

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

    node.deserialize(in);

1037
1038
1039

    // === Read admins. ===

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

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

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

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

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

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

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

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

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

	  index = elIndexVecIndex[index];

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

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

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

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

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

    initialized = true;
  }

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

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

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

  void Mesh::clearMacroFileInfo()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1185
    macroFileInfo->clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
1186
    delete macroFileInfo;
1187
    macroFileInfo = NULL;
1188
  }
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
  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;
    }
  }

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

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