Mesh.cc 33.1 KB
Newer Older
1001
    int oldVal = dim;
1002
1003
    SerUtil::deserialize(in, dim);
    TEST_EXIT_DBG(oldVal == 0 || dim == oldVal)("Invalid dimension!\n");
1004

1005
1006
1007
1008
1009
1010
    SerUtil::deserialize(in, nVertices);
    SerUtil::deserialize(in, nEdges);
    SerUtil::deserialize(in, nLeaves);
    SerUtil::deserialize(in, nElements);
    SerUtil::deserialize(in, nFaces);
    SerUtil::deserialize(in, maxEdgeNeigh);
1011
1012
1013

    diam.deserialize(in);

1014
    SerUtil::deserialize(in, preserveCoarseDOFs);
1015
1016

    oldVal = nDOFEl;
1017
1018
    SerUtil::deserialize(in, nDOFEl);
    TEST_EXIT_DBG(oldVal == 0 || nDOFEl == oldVal)("Invalid nDOFEl!\n");
1019
1020
1021
1022

    nDOF.deserialize(in);

    oldVal = nNodeEl;
1023
1024
    SerUtil::deserialize(in, nNodeEl);
    TEST_EXIT_DBG(oldVal == 0 || nNodeEl == oldVal)("Invalid nNodeEl!\n");
1025
1026
1027

    node.deserialize(in);

1028
1029
1030

    // === Read admins. ===

1031
    int size;
1032
    SerUtil::deserialize(in, size);
1033
    admin.resize(size, NULL);
1034
    for (int i = 0; i < size; i++) {
1035
      if (!admin[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1036
	admin[i] = new DOFAdmin(this);
1037

1038
1039
1040
      admin[i]->deserialize(in);
    }

1041
    SerUtil::deserialize(in, size);
Thomas Witkowski's avatar
Thomas Witkowski committed
1042
    std::vector< std::vector<int> > neighbourIndices(size);
1043

1044
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++)
1045
      if (macroElements[i])
Thomas Witkowski's avatar
Thomas Witkowski committed
1046
	delete macroElements[i];
1047
1048
1049
1050
1051
1052
1053
1054

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

1056
    macroElements.resize(size);
1057
    for (int i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1058
      macroElements[i] = new MacroElement(dim);
1059
1060
      macroElements[i]->writeNeighboursTo(&(neighbourIndices[i]));
      macroElements[i]->deserialize(in);
1061
      elIndexVecIndex[macroElements[i]->getIndex()] = i;
1062
1063
    }

1064
1065
    SerUtil::deserialize(in, elementIndex);
    SerUtil::deserialize(in, initialized);
1066
1067

    // set neighbour pointer in macro elements
1068
    int nNeighbour = getGeo(NEIGH);
1069
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++) {
1070
      for (int j = 0; j < nNeighbour; j++) {
1071
	int index = neighbourIndices[i][j];
1072
1073
1074
1075
1076
1077
1078

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

	  index = elIndexVecIndex[index];

1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
	  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);
1089
    while (elInfo) {
1090
1091
1092
1093
1094
      elInfo->getElement()->setMesh(this);
      elInfo = stack.traverseNext(elInfo);
    }

    serializedDOFs.clear();
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110

    
    /// === 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;      
    }
1111
1112
1113
1114
  }

  void Mesh::initialize() 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1115
1116
1117
    std::string macroFilename("");
    std::string valueFilename("");
    std::string periodicFile("");
1118
1119
1120
1121
1122
1123
1124
1125
1126
    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()) {
1127
      macroFileInfo = MacroReader::readMacro(macroFilename.c_str(), this,
1128
1129
 					     periodicFile == "" ? NULL : periodicFile.c_str(),
 					     check);
1130
1131
1132

      // If there is no value file which should be written, we can delete
      // the information of the macro file.
1133
1134
      if (!valueFilename.length())
  	clearMacroFileInfo();
1135
1136
1137
1138
1139
    }

    initialized = true;
  }

1140
1141
  bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1142
1143
    std::map<BoundaryType, VertexVector*>::iterator it;
    std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
1144
    for (it = periodicAssociations.begin(); it != end; ++it)
1145
1146
1147
1148
1149
      if ((*(it->second))[dof1] == dof2)
	return true;
    return false;
  }

1150
1151
  bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1152
1153
1154
    std::vector<DegreeOfFreedom> associatedToDOF1;
    std::map<BoundaryType, VertexVector*>::iterator it;
    std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
1155
1156
1157
    DegreeOfFreedom dof, assDOF;

    associatedToDOF1.push_back(dof1);
Thomas Witkowski's avatar
Thomas Witkowski committed
1158
1159
1160
    for (it = periodicAssociations.begin(); it != end; ++it) {
      int size = static_cast<int>(associatedToDOF1.size());
      for (int i = 0; i < size; i++) {
1161
1162
	dof = associatedToDOF1[i];
	assDOF = (*(it->second))[dof];
Thomas Witkowski's avatar
Thomas Witkowski committed
1163
	if (assDOF == dof2) {
1164
1165
	  return true;
	} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1166
1167
	  if (assDOF != dof) 
	    associatedToDOF1.push_back(assDOF);
1168
1169
1170
1171
1172
1173
1174
1175
	}
      }
    }
    return false;
  }

  void Mesh::clearMacroFileInfo()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1176
    macroFileInfo->clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
1177
    delete macroFileInfo;
1178
    macroFileInfo = NULL;
1179
  }
1180

1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
  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;
    }
  }

1209
1210
  int Mesh::calcMemoryUsage()
  {
1211
    int result = sizeof(Mesh);
1212

1213
1214
1215
1216
1217
1218
    result += nDOFEl;
    for (int i = 0; i < static_cast<int>(admin.size()); i++) {
      result += admin[i]->calcMemoryUsage();
      result += admin[i]->getUsedSize() * sizeof(DegreeOfFreedom);
    }
    
1219
1220
    for (int i = 0; i < static_cast<int>(macroElements.size()); i++)
      result += macroElements[i]->calcMemoryUsage();    
1221
1222
1223
    
    return result;
  }
1224
}
For faster browsing, not all history is shown. View entire blame