// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to <gmsh@geuz.org>. #include "PViewDataList.h" #include "GmshMessage.h" #include "GmshDefines.h" #include "polynomialBasis.h" #include "Numeric.h" #include "SmoothData.h" #include "Context.h" PViewDataList::PViewDataList(bool isAdapted) : PViewData(), NbTimeStep(0), Min(VAL_INF), Max(-VAL_INF), NbSP(0), NbVP(0), NbTP(0), NbSL(0), NbVL(0), NbTL(0), NbST(0), NbVT(0), NbTT(0), NbSQ(0), NbVQ(0), NbTQ(0), NbSG(0), NbVG(0), NbTG(0), NbSS(0), NbVS(0), NbTS(0), NbSH(0), NbVH(0), NbTH(0), NbSI(0), NbVI(0), NbTI(0), NbSY(0), NbVY(0), NbTY(0), NbSD(0), NbVD(0), NbTD(0), NbT2(0), NbT3(0), _lastElement(-1), _lastDimension(-1), _lastNumNodes(-1), _lastNumComponents(-1), _lastNumValues(-1), _lastNumEdges(-1), _lastType(-1), _lastXYZ(0), _lastVal(0), _isAdapted(isAdapted) { for(int i = 0; i < 30; i++) _index[i] = 0; polyTotNumNodes[0] = 0.; polyTotNumNodes[1] = 0.; polyAgNumNodes[0].push_back(0.); polyAgNumNodes[1].push_back(0.); } void PViewDataList::setXY(std::vector<double> &x, std::vector<double> &y) { NbSP = 0; SP.clear(); for(unsigned int i = 0; i < std::min(x.size(), y.size()); i++){ SP.push_back(x[i]); SP.push_back(0.); SP.push_back(0.); SP.push_back(y[i]); NbSP++; } finalize(); } bool PViewDataList::finalize(bool computeMinMax, const std::string &interpolationScheme) { BBox.reset(); Min = VAL_INF; Max = -VAL_INF; // finalize text strings first, to get the max value of NbTimeStep // for strings-only views (strings are designed to degrade // gracefully when some have fewer time steps than others). If there // are any elements in the view, this value will be replaced by the // minimum number of time steps common to all elements. _stat(T2D, T2C, 4); _stat(T3D, T3C, 5); // compute min/max and other statistics for all element lists _stat(SP, 1, NbSP, 1, TYPE_PNT); _stat(VP, 3, NbVP, 1, TYPE_PNT); _stat(TP, 9, NbTP, 1, TYPE_PNT); _stat(SL, 1, NbSL, 2, TYPE_LIN); _stat(VL, 3, NbVL, 2, TYPE_LIN); _stat(TL, 9, NbTL, 2, TYPE_LIN); _stat(ST, 1, NbST, 3, TYPE_TRI); _stat(VT, 3, NbVT, 3, TYPE_TRI); _stat(TT, 9, NbTT, 3, TYPE_TRI); _stat(SQ, 1, NbSQ, 4, TYPE_QUA); _stat(VQ, 3, NbVQ, 4, TYPE_QUA); _stat(TQ, 9, NbTQ, 4, TYPE_QUA); _stat(SS, 1, NbSS, 4, TYPE_TET); _stat(VS, 3, NbVS, 4, TYPE_TET); _stat(TS, 9, NbTS, 4, TYPE_TET); _stat(SH, 1, NbSH, 8, TYPE_HEX); _stat(VH, 3, NbVH, 8, TYPE_HEX); _stat(TH, 9, NbTH, 8, TYPE_HEX); _stat(SI, 1, NbSI, 6, TYPE_PRI); _stat(VI, 3, NbVI, 6, TYPE_PRI); _stat(TI, 9, NbTI, 6, TYPE_PRI); _stat(SY, 1, NbSY, 5, TYPE_PYR); _stat(VY, 3, NbVY, 5, TYPE_PYR); _stat(TY, 9, NbTY, 5, TYPE_PYR); _stat(SG, 1, NbSG, 3, TYPE_POLYG); _stat(VG, 3, NbVG, 3, TYPE_POLYG); _stat(TG, 9, NbTG, 3, TYPE_POLYG); _stat(SD, 1, NbSD, 4, TYPE_POLYH); _stat(VD, 3, NbVD, 4, TYPE_POLYH); _stat(TD, 9, NbTD, 4, TYPE_POLYH); // add dummy time values if none (or too few) time values are // provided (e.g. using the old parsed format) if((int)Time.size() < NbTimeStep) { for(int i = Time.size(); i < NbTimeStep; i++) Time.push_back(i); } // compute starting element indices int nb[30] = {NbSP, NbVP, NbTP, NbSL, NbVL, NbTL, NbST, NbVT, NbTT, NbSQ, NbVQ, NbTQ, NbSS, NbVS, NbTS, NbSH, NbVH, NbTH, NbSI, NbVI, NbTI, NbSY, NbVY, NbTY, NbSG, NbVG, NbTG, NbSD, NbVD, NbTD}; for(int i = 0; i < 30; i++){ _index[i] = 0; for(int j = 0; j <= i; j++) _index[i] += nb[j]; } if(CTX::instance()->post.smooth) smooth(); return PViewData::finalize(); } int PViewDataList::getNumScalars(int step) { return NbSP + NbSL + NbST + NbSQ + NbSS + NbSH + NbSI + NbSY + NbSG + NbSD; } int PViewDataList::getNumVectors(int step) { return NbVP + NbVL + NbVT + NbVQ + NbVS + NbVH + NbVI + NbVY + NbVG + NbVD; } int PViewDataList::getNumTensors(int step) { return NbTP + NbTL + NbTT + NbTQ + NbTS + NbTH + NbTI + NbTY + NbTG + NbTD; } int PViewDataList::getNumElements(int step, int ent) { return getNumScalars() + getNumVectors() + getNumTensors(); } double PViewDataList::getTime(int step) { if(step < 0 || step >= (int)Time.size()) return 0.; return Time[step]; } double PViewDataList::getMin(int step, bool onlyVisible, int forceNumComponents, int componentMap[9]) { if(step >= (int)TimeStepMin.size()) return Min; if(forceNumComponents){ double vmin = VAL_INF; for(int ent = 0; ent < getNumEntities(step); ent++){ for(int ele = 0; ele < getNumElements(step, ent); ele++){ for(int nod = 0; nod < getNumNodes(step, ent, ele); nod++){ double val; getScalarValue(step, ent, ele, nod, val, forceNumComponents, componentMap); vmin = std::min(vmin, val); } } } return vmin; } if(step < 0) return Min; return TimeStepMin[step]; } double PViewDataList::getMax(int step, bool onlyVisible, int forceNumComponents, int componentMap[9]) { if(step >= (int)TimeStepMax.size()) return Max; if(forceNumComponents){ double vmax = -VAL_INF; for(int ent = 0; ent < getNumEntities(step); ent++){ for(int ele = 0; ele < getNumElements(step, ent); ele++){ for(int nod = 0; nod < getNumNodes(step, ent, ele); nod++){ double val; getScalarValue(step, ent, ele, nod, val, forceNumComponents, componentMap); vmax = std::max(vmax, val); } } } return vmax; } if(step < 0) return Max; return TimeStepMax[step]; } void PViewDataList::_stat(std::vector<double> &D, std::vector<char> &C, int nb) { // compute statistics for text lists for(unsigned int i = 0; i < D.size(); i += nb){ double beg = D[i + nb - 1]; double end; if(i + 2 * nb > D.size()) end = C.size(); else end = D[i + nb + nb - 1]; char *c = &C[(int)beg]; int nbtime = 0; for(int j = 0; j < (int)(end - beg); j++) if(c[j] == '\0') nbtime++; if(nbtime > NbTimeStep) NbTimeStep = nbtime; } if(nb == 5){ for(unsigned int i = 0; i < D.size(); i += nb) BBox += SPoint3(D[i], D[i + 1], D[i + 2]); } } void PViewDataList::_stat(std::vector<double> &list, int nbcomp, int nbelm, int nbnod, int type) { // compute statistics for element lists if(!nbelm) return; int nbval = nbcomp * nbnod; if(haveInterpolationMatrices()){ std::vector<fullMatrix<double>*> im; int nim = getInterpolationMatrices(type, im); if(nim == 4) nbnod = im[2]->size1(); if(nim) nbval = nbcomp * im[0]->size1(); } int nb = list.size() / nbelm; for(int ele = 0; ele < nbelm; ele ++){ int i = ele * nb; if(type == TYPE_POLYG || type == TYPE_POLYH){ int t = (type == TYPE_POLYG) ? 0 : 1; nbnod = polyNumNodes[t][ele]; nb = list.size() / polyTotNumNodes[t] * nbnod; i = polyAgNumNodes[t][ele] * nb / nbnod; nbval = nbcomp * nbnod; } int N = nb - 3 * nbnod; double *X = &list[i]; double *Y = &list[i + 1 * nbnod]; double *Z = &list[i + 2 * nbnod]; double *V = &list[i + 3 * nbnod]; // update bounding box for(int j = 0; j < nbnod; j++) BBox += SPoint3(X[j], Y[j], Z[j]); // update num time steps if(Min == VAL_INF || Max == -VAL_INF){ NbTimeStep = N / nbval; TimeStepMin.clear(); TimeStepMax.clear(); for(int j = 0; j < NbTimeStep; j++){ TimeStepMin.push_back(VAL_INF); TimeStepMax.push_back(-VAL_INF); } } else if(N / nbval < NbTimeStep){ // if some elts have less steps, reduce the total number! NbTimeStep = N / nbval; } // update min/max for(int j = 0; j < N; j += nbcomp) { double l0 = ComputeScalarRep(nbcomp, &V[j]); Min = std::min(l0, Min); Max = std::max(l0, Max); int ts = j / nbval; if(ts < NbTimeStep){ // security TimeStepMin[ts] = std::min(l0, TimeStepMin[ts]); TimeStepMax[ts] = std::max(l0, TimeStepMax[ts]); } } } } void PViewDataList::_setLast(int ele, int dim, int nbnod, int nbcomp, int nbedg, int type, std::vector<double> &list, int nblist) { if(haveInterpolationMatrices()){ std::vector<fullMatrix<double>*> im; if(getInterpolationMatrices(type, im) == 4) nbnod = im[2]->size1(); } _lastDimension = dim; _lastNumNodes = nbnod; _lastNumComponents = nbcomp; _lastNumEdges = nbedg; _lastType = type; int nb = list.size() / nblist; // number of coords and values for the element int nbAg = ele * nb; // number of coords and values before the ones of the element if(type == TYPE_POLYG || type == TYPE_POLYH){ int t = (type == TYPE_POLYG) ? 0 : 1; nb = list.size() / polyTotNumNodes[t] * nbnod; nbAg = polyAgNumNodes[t][ele] * nb / nbnod; } _lastNumValues = (nb - 3 * nbnod) / NbTimeStep; _lastXYZ = &list[nbAg]; _lastVal = &list[nbAg + 3 * _lastNumNodes]; } void PViewDataList::_setLast(int ele) { _lastElement = ele; if(ele < _index[2]){ // points if(ele < _index[0]) _setLast(ele, 0, 1, 1, 0, TYPE_PNT, SP, NbSP); else if(ele < _index[1]) _setLast(ele - _index[0], 0, 1, 3, 0, TYPE_PNT, VP, NbVP); else _setLast(ele - _index[1], 0, 1, 9, 0, TYPE_PNT, TP, NbTP); } else if(ele < _index[5]){ // lines if(ele < _index[3]) _setLast(ele - _index[2], 1, 2, 1, 1, TYPE_LIN, SL, NbSL); else if(ele < _index[4]) _setLast(ele - _index[3], 1, 2, 3, 1, TYPE_LIN, VL, NbVL); else _setLast(ele - _index[4], 1, 2, 9, 1, TYPE_LIN, TL, NbTL); } else if(ele < _index[8]){ // triangles if(ele < _index[6]) _setLast(ele - _index[5], 2, 3, 1, 3, TYPE_TRI, ST, NbST); else if(ele < _index[7]) _setLast(ele - _index[6], 2, 3, 3, 3, TYPE_TRI, VT, NbVT); else _setLast(ele - _index[7], 2, 3, 9, 3, TYPE_TRI, TT, NbTT); } else if(ele < _index[11]){ // quadrangles if(ele < _index[9]) _setLast(ele - _index[8], 2, 4, 1, 4, TYPE_QUA, SQ, NbSQ); else if(ele < _index[10]) _setLast(ele - _index[9], 2, 4, 3, 4, TYPE_QUA, VQ, NbVQ); else _setLast(ele - _index[10], 2, 4, 9, 4, TYPE_QUA, TQ, NbTQ); } else if(ele < _index[14]){ // tetrahedra if(ele < _index[12]) _setLast(ele - _index[11], 3, 4, 1, 6, TYPE_TET, SS, NbSS); else if(ele < _index[13]) _setLast(ele - _index[12], 3, 4, 3, 6, TYPE_TET, VS, NbVS); else _setLast(ele - _index[13], 3, 4, 9, 6, TYPE_TET, TS, NbTS); } else if(ele < _index[17]){ // hexahedra if(ele < _index[15]) _setLast(ele - _index[14], 3, 8, 1, 12, TYPE_HEX, SH, NbSH); else if(ele < _index[16]) _setLast(ele - _index[15], 3, 8, 3, 12, TYPE_HEX, VH, NbVH); else _setLast(ele - _index[16], 3, 8, 9, 12, TYPE_HEX, TH, NbTH); } else if(ele < _index[20]){ // prisms if(ele < _index[18]) _setLast(ele - _index[17], 3, 6, 1, 9, TYPE_PRI, SI, NbSI); else if(ele < _index[19]) _setLast(ele - _index[18], 3, 6, 3, 9, TYPE_PRI, VI, NbVI); else _setLast(ele - _index[19], 3, 6, 9, 9, TYPE_PRI, TI, NbTI); } else if(ele < _index[23]){ // pyramids if(ele < _index[21]) _setLast(ele - _index[20], 3, 5, 1, 8, TYPE_PYR, SY, NbSY); else if(ele < _index[22]) _setLast(ele - _index[21], 3, 5, 3, 8, TYPE_PYR, VY, NbVY); else _setLast(ele - _index[22], 3, 5, 9, 8, TYPE_PYR, TY, NbTY); } else if(ele < _index[26]){ // polygons int nN = polyNumNodes[0][ele - _index[23]]; if(ele < _index[24]) _setLast(ele - _index[23], 2, nN, 1, nN, TYPE_POLYG, SG, NbSG); else if(ele < _index[25]) _setLast(ele - _index[24], 2, nN, 3, nN, TYPE_POLYG, VG, NbVG); else _setLast(ele - _index[25], 2, nN, 9, nN, TYPE_POLYG, TG, NbTG); } else if(ele < _index[29]){ // polyhedra int nN = polyNumNodes[1][ele - _index[26]]; if(ele < _index[27]) _setLast(ele - _index[26], 3, nN, 1, nN*1.5, TYPE_POLYH, SD, NbSD); else if(ele < _index[28]) _setLast(ele - _index[27], 3, nN, 3, nN*1.5, TYPE_POLYH, VD, NbVD); else _setLast(ele - _index[28], 3, nN, 9, nN*1.5, TYPE_POLYH, TD, NbTD); } } int PViewDataList::getDimension(int step, int ent, int ele) { if(ele != _lastElement) _setLast(ele); return _lastDimension; } int PViewDataList::getNumNodes(int step, int ent, int ele) { if(ele != _lastElement) _setLast(ele); return _lastNumNodes; } int PViewDataList::getNode(int step, int ent, int ele, int nod, double &x, double &y, double &z) { if(ele != _lastElement) _setLast(ele); x = _lastXYZ[nod]; y = _lastXYZ[_lastNumNodes + nod]; z = _lastXYZ[2 * _lastNumNodes + nod]; return 0; } void PViewDataList::setNode(int step, int ent, int ele, int nod, double x, double y, double z) { if(step) return; if(ele != _lastElement) _setLast(ele); _lastXYZ[nod] = x; _lastXYZ[_lastNumNodes + nod] = y; _lastXYZ[2 * _lastNumNodes + nod] = z; } int PViewDataList::getNumComponents(int step, int ent, int ele) { if(ele != _lastElement) _setLast(ele); return _lastNumComponents; } int PViewDataList::getNumValues(int step, int ent, int ele) { if(ele != _lastElement) _setLast(ele); return _lastNumValues; } void PViewDataList::getValue(int step, int ent, int ele, int idx, double &val) { if(ele != _lastElement) _setLast(ele); if(step >= NbTimeStep) step = 0; val = _lastVal[step * _lastNumValues + idx]; } void PViewDataList::getValue(int step, int ent, int ele, int nod, int comp, double &val) { if(ele != _lastElement) _setLast(ele); if(step >= NbTimeStep) step = 0; val = _lastVal[step * _lastNumNodes * _lastNumComponents + nod * _lastNumComponents + comp]; } void PViewDataList::setValue(int step, int ent, int ele, int nod, int comp, double val) { if(ele != _lastElement) _setLast(ele); if(step >= NbTimeStep) step = 0; _lastVal[step * _lastNumNodes * _lastNumComponents + nod * _lastNumComponents + comp] = val; } int PViewDataList::getNumEdges(int step, int ent, int ele) { if(ele != _lastElement) _setLast(ele); return _lastNumEdges; } int PViewDataList::getType(int step, int ent, int ele) { if(ele != _lastElement) _setLast(ele); return _lastType; } void PViewDataList::_getString(int dim, int i, int step, std::string &str, double &x, double &y, double &z, double &style) { // 3D: T3D is a list of double: x,y,z,style,index,x,y,z,style,index,... // T3C is a list of chars: string\0,string\0,string\0,string\0,... // Parser format is: T3(x,y,z,style){"str","str",...}; // 2D: T2D is a list of double: x,y,style,index,x,y,style,index,... // T2C is a list of chars: string\0,string\0,string\0,string\0,... // Parser format is: T2(x,y,style){"str","str",...}; std::vector<double> &td = (dim == 2) ? T2D : T3D; std::vector<char> &tc = (dim == 2) ? T2C : T3C; int nbd = (dim == 2) ? 4 : 5; int index, nbchar; double *d1 = &td[i * nbd]; double *d2 = ((i + 1) * nbd < (int)td.size()) ? &td[(i + 1) * nbd] : 0; if(dim == 2) { x = d1[0]; y = d1[1]; z = 0.; style = d1[2]; index = (int)d1[3]; if(d2) nbchar = (int)d2[3] - index; else nbchar = tc.size() - index; } else { x = d1[0]; y = d1[1]; z = d1[2]; style = d1[3]; index = (int)d1[4]; if(d2) nbchar = (int)d2[4] - index; else nbchar = tc.size() - index; } char *c = &tc[index]; int k = 0, l = 0; while(k < nbchar && l != step) { if(c[k++] == '\0') l++; } if(k < nbchar && l == step) str = std::string(&c[k]); else str = std::string(c); } void PViewDataList::getString2D(int i, int step, std::string &str, double &x, double &y, double &style) { double z; _getString(2, i, step, str, x, y, z, style); } void PViewDataList::getString3D(int i, int step, std::string &str, double &x, double &y, double &z, double &style) { _getString(3, i, step, str, x, y, z, style); } void PViewDataList::revertElement(int step, int ent, int ele) { if(step) return; if(ele != _lastElement) _setLast(ele); // copy data std::vector<double> XYZ(3 * _lastNumNodes); for(unsigned int i = 0; i < XYZ.size(); i++) XYZ[i] = _lastXYZ[i]; std::vector<double> V(_lastNumNodes * _lastNumComponents * getNumTimeSteps()); for(unsigned int i = 0; i < V.size(); i++) V[i] = _lastVal[i]; // reverse node order for(int i = 0; i < _lastNumNodes; i++){ _lastXYZ[i] = XYZ[_lastNumNodes - i - 1]; _lastXYZ[_lastNumNodes + i] = XYZ[2 * _lastNumNodes - i - 1]; _lastXYZ[2 * _lastNumNodes + i] = XYZ[3 * _lastNumNodes - i - 1]; } for(int step = 0; step < getNumTimeSteps(); step++) for(int i = 0; i < _lastNumNodes; i++) for(int k = 0; k < _lastNumComponents; k++) _lastVal[_lastNumComponents * _lastNumNodes * step + _lastNumComponents * i + k] = V[_lastNumComponents * _lastNumNodes * step + _lastNumComponents * (_lastNumNodes - i - 1) + k]; } static void generateConnectivities(std::vector<double> &list, int nbList, int nbTimeStep, int nbVert, int nbComp, smooth_data &data) { if(!nbList) return; double *vals = new double[nbTimeStep * nbComp]; int nb = list.size() / nbList; for(unsigned int i = 0; i < list.size(); i += nb) { double *x = &list[i]; double *y = &list[i + nbVert]; double *z = &list[i + 2 * nbVert]; double *v = &list[i + 3 * nbVert]; for(int j = 0; j < nbVert; j++) { for(int ts = 0; ts < nbTimeStep; ts++) for(int k = 0; k < nbComp; k++) vals[nbComp * ts + k] = v[nbVert * nbComp * ts + nbComp * j + k]; data.add(x[j], y[j], z[j], nbTimeStep * nbComp, vals); } } delete [] vals; } static void smoothList(std::vector<double> &list, int nbList, int nbTimeStep, int nbVert, int nbComp, smooth_data &data) { if(!nbList) return; double *vals = new double[nbTimeStep * nbComp]; int nb = list.size() / nbList; for(unsigned int i = 0; i < list.size(); i += nb) { double *x = &list[i]; double *y = &list[i + nbVert]; double *z = &list[i + 2 * nbVert]; double *v = &list[i + 3 * nbVert]; for(int j = 0; j < nbVert; j++) { if(data.get(x[j], y[j], z[j], nbTimeStep * nbComp, vals)){ for(int ts = 0; ts < nbTimeStep; ts++) for(int k = 0; k < nbComp; k++) v[nbVert * nbComp * ts + nbComp * j + k] = vals[nbComp * ts + k]; } } } delete [] vals; } void PViewDataList::smooth() { double old_eps = xyzv::eps; xyzv::eps = CTX::instance()->lc * 1.e-8; smooth_data data; std::vector<double> *list = 0; int *nbe = 0, nbc, nbn; for(int i = 0; i < 24; i++){ _getRawData(i, &list, &nbe, &nbc, &nbn); if(nbn > 1) generateConnectivities(*list, *nbe, NbTimeStep, nbn, nbc, data); } for(int i = 0; i < 24; i++){ _getRawData(i, &list, &nbe, &nbc, &nbn); if(nbn > 1) smoothList(*list, *nbe, NbTimeStep, nbn, nbc, data); } xyzv::eps = old_eps; finalize(); } static void dVecMerge(std::vector<double> &v, std::vector<double> &dest) { for(unsigned int i = 0; i < v.size(); i++) dest.push_back(v[i]); } bool PViewDataList::combineSpace(nameData &nd) { // sanity checks if(nd.data.size() < 2) return false; int ts = nd.data[0]->getNumTimeSteps(); for(unsigned int i = 1; i < nd.data.size(); i++) { if(!nd.data[i]->empty() && nd.data[i]->getNumTimeSteps() != ts){ Msg::Error("Cannot combine views having different number of time steps"); return false; } } for(unsigned int i = 0; i < nd.data.size(); i++) { PViewDataList *l = dynamic_cast<PViewDataList*>(nd.data[i]); if(!l){ Msg::Error("Cannot combine hybrid data"); return false; } // copy interpolation marices for(std::map<int, std::vector<fullMatrix<double>*> >::iterator it = l->_interpolation.begin(); it != l->_interpolation.end(); it++) if(_interpolation[it->first].empty()) for(unsigned int i = 0; i < it->second.size(); i++) _interpolation[it->first].push_back(new fullMatrix<double>(*it->second[i])); // merge elememts dVecMerge(l->SP, SP); NbSP += l->NbSP; dVecMerge(l->VP, VP); NbVP += l->NbVP; dVecMerge(l->TP, TP); NbTP += l->NbTP; dVecMerge(l->SL, SL); NbSL += l->NbSL; dVecMerge(l->VL, VL); NbVL += l->NbVL; dVecMerge(l->TL, TL); NbTL += l->NbTL; dVecMerge(l->ST, ST); NbST += l->NbST; dVecMerge(l->VT, VT); NbVT += l->NbVT; dVecMerge(l->TT, TT); NbTT += l->NbTT; dVecMerge(l->SQ, SQ); NbSQ += l->NbSQ; dVecMerge(l->VQ, VQ); NbVQ += l->NbVQ; dVecMerge(l->TQ, TQ); NbTQ += l->NbTQ; dVecMerge(l->SS, SS); NbSS += l->NbSS; dVecMerge(l->VS, VS); NbVS += l->NbVS; dVecMerge(l->TS, TS); NbTS += l->NbTS; dVecMerge(l->SH, SH); NbSH += l->NbSH; dVecMerge(l->VH, VH); NbVH += l->NbVH; dVecMerge(l->TH, TH); NbTH += l->NbTH; dVecMerge(l->SI, SI); NbSI += l->NbSI; dVecMerge(l->VI, VI); NbVI += l->NbVI; dVecMerge(l->TI, TI); NbTI += l->NbTI; dVecMerge(l->SY, SY); NbSY += l->NbSY; dVecMerge(l->VY, VY); NbVY += l->NbVY; dVecMerge(l->TY, TY); NbTY += l->NbTY; // merge strings for(unsigned int i = 0; i < l->T2D.size(); i += 4){ T2D.push_back(l->T2D[i]); T2D.push_back(l->T2D[i + 1]); T2D.push_back(l->T2D[i + 2]); T2D.push_back(T2C.size()); double beg = l->T2D[i + 3]; double end; if(i > l->T2D.size() - 8) end = l->T2C.size(); else end = l->T2D[i + 3 + 4]; char *c = &l->T2C[(int)beg]; for(int j = 0; j < (int)(end - beg); j++) T2C.push_back(c[j]); NbT2++; } for(unsigned int i = 0; i < l->T3D.size(); i += 5){ T3D.push_back(l->T3D[i]); T3D.push_back(l->T3D[i + 1]); T3D.push_back(l->T3D[i + 2]); T3D.push_back(l->T3D[i + 3]); T3D.push_back(T3C.size()); double beg = l->T3D[i + 4]; double end; if(i > l->T3D.size() - 10) end = l->T3C.size(); else end = l->T3D[i + 4 + 5]; char *c = &l->T3C[(int)beg]; for(int j = 0; j < (int)(end-beg); j++) T3C.push_back(c[j]); NbT3++; } } std::string tmp; if(nd.name == "__all__") tmp = "all"; else if(nd.name == "__vis__") tmp = "visible"; else tmp = nd.name; char name[256]; sprintf(name, "%s_Combine", tmp.c_str()); setName(name); setFileName(std::string(name) + ".pos"); return finalize(); } bool PViewDataList::combineTime(nameData &nd) { // sanity checks if(nd.data.size() < 2) return false; std::vector<PViewDataList*> data(nd.data.size()); for(unsigned int i = 0; i < nd.data.size(); i++){ data[i] = dynamic_cast<PViewDataList*>(nd.data[i]); if(!data[i]){ Msg::Error("Cannot combine hybrid data"); return false; } } int *nbe = 0, *nbe2 = 0, nbn, nbn2, nbc, nbc2; std::vector<double> *list = 0, *list2 = 0; // use the first data set as the reference for(int i = 0; i < 24; i++){ _getRawData(i, &list, &nbe, &nbc, &nbn); data[0]->_getRawData(i, &list2, &nbe2, &nbc2, &nbn2); *nbe = *nbe2; } NbT2 = data[0]->NbT2; NbT3 = data[0]->NbT3; for(std::map<int, std::vector<fullMatrix<double>*> >::iterator it = data[0]->_interpolation.begin(); it != data[0]->_interpolation.end(); it++) if(_interpolation[it->first].empty()) for(unsigned int i = 0; i < it->second.size(); i++) _interpolation[it->first].push_back(new fullMatrix<double>(*it->second[i])); // merge values for all element types for(int i = 0; i < 24; i++){ _getRawData(i, &list, &nbe, &nbc, &nbn); for(int j = 0; j < *nbe; j++){ for(unsigned int k = 0; k < data.size(); k++){ data[k]->_getRawData(i, &list2, &nbe2, &nbc2, &nbn2); if(*nbe && *nbe == *nbe2){ int nb2 = list2->size() / *nbe2; if(!k){ // copy coordinates of elm j (we are always here as // expected, since the ref view is the first one) for(int l = 0; l < 3 * nbn2; l++) list->push_back((*list2)[j * nb2 + l]); } // copy values of elm j for(int l = 0; l < nb2 - 3 * nbn2; l++) list->push_back((*list2)[j * nb2 + 3 * nbn2 + l]); } } } } // merge 2d strings for(int j = 0; j < NbT2; j++){ for(unsigned int k = 0; k < data.size(); k++){ if(NbT2 == data[k]->NbT2){ if(!k){ // copy coordinates T2D.push_back(data[k]->T2D[j * 4]); T2D.push_back(data[k]->T2D[j * 4 + 1]); T2D.push_back(data[k]->T2D[j * 4 + 2]); // index T2D.push_back(T2C.size()); } // copy char values double beg = data[k]->T2D[j * 4 + 3]; double end; if(j == NbT2 - 1) end = data[k]->T2C.size(); else end = data[k]->T2D[j * 4 + 4 + 3]; char *c = &data[k]->T2C[(int)beg]; for(int l = 0; l < (int)(end - beg); l++) T2C.push_back(c[l]); } } } // merge 3d strings for(int j = 0; j < NbT3; j++){ for(unsigned int k = 0; k < data.size(); k++){ if(NbT3 == data[k]->NbT3){ if(!k){ // copy coordinates T3D.push_back(data[k]->T3D[j * 5]); T3D.push_back(data[k]->T3D[j * 5 + 1]); T3D.push_back(data[k]->T3D[j * 5 + 2]); T3D.push_back(data[k]->T3D[j * 5 + 3]); // index T3D.push_back(T3C.size()); } // copy char values double beg = data[k]->T3D[j * 5 + 4]; double end; if(j == NbT3 - 1) end = data[k]->T3C.size(); else end = data[k]->T3D[j * 5 + 5 + 4]; char *c = &data[k]->T3C[(int)beg]; for(int l = 0; l < (int)(end - beg); l++) T3C.push_back(c[l]); } } } // create the time data for(unsigned int i = 0; i < data.size(); i++) dVecMerge(data[i]->Time, Time); // if all the time values are the same, it probably means that the // original views didn't have any time data: then we'll just use // time step values if(Time.size()){ double t0 = Time[0], ti; bool allTheSame = true; for(unsigned int i = 1; i < Time.size(); i++){ ti = Time[i]; if(ti != t0){ allTheSame = false; break; } } if(allTheSame) Time.clear(); } std::string tmp; if(nd.name == "__all__") tmp = "all"; else if(nd.name == "__vis__") tmp = "visible"; else tmp = nd.name; char name[256]; sprintf(name, "%s_Combine", tmp.c_str()); setName(name); setFileName(std::string(name) + ".pos"); return finalize(); } int PViewDataList::_getRawData(int idxtype, std::vector<double> **l, int **ne, int *nc, int *nn) { int type = 0; // No constant nn for polygons! if(idxtype > 23 && idxtype < 30) Msg::Warning("No constant number of nodes for polygons and polyhedra."); switch(idxtype){ case 0 : *l = &SP; *ne = &NbSP; *nc = 1; *nn = 1; type = TYPE_PNT; break; case 1 : *l = &VP; *ne = &NbVP; *nc = 3; *nn = 1; type = TYPE_PNT; break; case 2 : *l = &TP; *ne = &NbTP; *nc = 9; *nn = 1; type = TYPE_PNT; break; case 3 : *l = &SL; *ne = &NbSL; *nc = 1; *nn = 2; type = TYPE_LIN; break; case 4 : *l = &VL; *ne = &NbVL; *nc = 3; *nn = 2; type = TYPE_LIN; break; case 5 : *l = &TL; *ne = &NbTL; *nc = 9; *nn = 2; type = TYPE_LIN; break; case 6 : *l = &ST; *ne = &NbST; *nc = 1; *nn = 3; type = TYPE_TRI; break; case 7 : *l = &VT; *ne = &NbVT; *nc = 3; *nn = 3; type = TYPE_TRI; break; case 8 : *l = &TT; *ne = &NbTT; *nc = 9; *nn = 3; type = TYPE_TRI; break; case 9 : *l = &SQ; *ne = &NbSQ; *nc = 1; *nn = 4; type = TYPE_QUA; break; case 10: *l = &VQ; *ne = &NbVQ; *nc = 3; *nn = 4; type = TYPE_QUA; break; case 11: *l = &TQ; *ne = &NbTQ; *nc = 9; *nn = 4; type = TYPE_QUA; break; case 12: *l = &SS; *ne = &NbSS; *nc = 1; *nn = 4; type = TYPE_TET; break; case 13: *l = &VS; *ne = &NbVS; *nc = 3; *nn = 4; type = TYPE_TET; break; case 14: *l = &TS; *ne = &NbTS; *nc = 9; *nn = 4; type = TYPE_TET; break; case 15: *l = &SH; *ne = &NbSH; *nc = 1; *nn = 8; type = TYPE_HEX; break; case 16: *l = &VH; *ne = &NbVH; *nc = 3; *nn = 8; type = TYPE_HEX; break; case 17: *l = &TH; *ne = &NbTH; *nc = 9; *nn = 8; type = TYPE_HEX; break; case 18: *l = &SI; *ne = &NbSI; *nc = 1; *nn = 6; type = TYPE_PRI; break; case 19: *l = &VI; *ne = &NbVI; *nc = 3; *nn = 6; type = TYPE_PRI; break; case 20: *l = &TI; *ne = &NbTI; *nc = 9; *nn = 6; type = TYPE_PRI; break; case 21: *l = &SY; *ne = &NbSY; *nc = 1; *nn = 5; type = TYPE_PYR; break; case 22: *l = &VY; *ne = &NbVY; *nc = 3; *nn = 5; type = TYPE_PYR; break; case 23: *l = &TY; *ne = &NbTY; *nc = 9; *nn = 5; type = TYPE_PYR; break; case 24: *l = &SG; *ne = &NbSG; *nc = 1; *nn = 3; type = TYPE_POLYG; break; case 25: *l = &VG; *ne = &NbVG; *nc = 3; *nn = 3; type = TYPE_POLYG; break; case 26: *l = &TG; *ne = &NbTG; *nc = 9; *nn = 3; type = TYPE_POLYG; break; case 27: *l = &SD; *ne = &NbSD; *nc = 1; *nn = 4; type = TYPE_POLYH; break; case 28: *l = &VD; *ne = &NbVD; *nc = 3; *nn = 4; type = TYPE_POLYH; break; case 29: *l = &TD; *ne = &NbTD; *nc = 9; *nn = 4; type = TYPE_POLYH; break; default: Msg::Error("Wrong type in PViewDataList"); break; } if(haveInterpolationMatrices()){ std::vector<fullMatrix<double>*> im; int nim = getInterpolationMatrices(type, im); if(nim == 4) *nn = im[2]->size1(); } return type; } void PViewDataList::setOrder2(int type) { int typeMSH = 0; switch(type){ case TYPE_LIN: typeMSH = MSH_LIN_3; break; case TYPE_TRI: typeMSH = MSH_TRI_6; break; case TYPE_QUA: typeMSH = MSH_QUA_9; break; case TYPE_TET: typeMSH = MSH_TET_10; break; case TYPE_HEX: typeMSH = MSH_HEX_27; break; case TYPE_PRI: typeMSH = MSH_PRI_18; break; case TYPE_PYR: typeMSH = MSH_PYR_14; break; } const polynomialBasis *fs = polynomialBases::find(typeMSH); if(!fs){ Msg::Error("Could not find function space for element type %d", typeMSH); return; } setInterpolationMatrices(type, fs->coefficients, fs->monomials, fs->coefficients, fs->monomials); } std::vector<double> *PViewDataList::incrementList(int numComp, int type, int numNodes) { int nb; switch(type){ case TYPE_PNT: if (numComp == 1){ NbSP++; return &SP; } else if(numComp == 3){ NbVP++; return &VP; } else if(numComp == 9){ NbTP++; return &TP; } break; case TYPE_LIN: if (numComp == 1){ NbSL++; return &SL; } else if(numComp == 3){ NbVL++; return &VL; } else if(numComp == 9){ NbTL++; return &TL; } break; case TYPE_TRI: if (numComp == 1){ NbST++; return &ST; } else if(numComp == 3){ NbVT++; return &VT; } else if(numComp == 9){ NbTT++; return &TT; } break; case TYPE_QUA: if (numComp == 1){ NbSQ++; return &SQ; } else if(numComp == 3){ NbVQ++; return &VQ; } else if(numComp == 9){ NbTQ++; return &TQ; } break; case TYPE_TET: if (numComp == 1){ NbSS++; return &SS; } else if(numComp == 3){ NbVS++; return &VS; } else if(numComp == 9){ NbTS++; return &TS; } break; case TYPE_HEX: if (numComp == 1){ NbSH++; return &SH; } else if(numComp == 3){ NbVH++; return &VH; } else if(numComp == 9){ NbTH++; return &TH; } break; case TYPE_PRI: if (numComp == 1){ NbSI++; return &SI; } else if(numComp == 3){ NbVI++; return &VI; } else if(numComp == 9){ NbTI++; return &TI; } break; case TYPE_PYR: if (numComp == 1){ NbSY++; return &SY; } else if(numComp == 3){ NbVY++; return &VY; } else if(numComp == 9){ NbTY++; return &TY; } break; case TYPE_POLYG: polyNumNodes[0].push_back(numNodes); nb = (polyAgNumNodes[0].size()) ? polyAgNumNodes[0].back() : 0; polyAgNumNodes[0].push_back(numNodes + nb); polyTotNumNodes[0] += numNodes; if (numComp == 1){ NbSG++; return &SG; } else if(numComp == 3){ NbVG++; return &VG; } else if(numComp == 9){ NbTG++; return &TG; } break; case TYPE_POLYH: polyNumNodes[1].push_back(numNodes); nb = (polyAgNumNodes[1].size()) ? polyAgNumNodes[1].back() : 0; polyAgNumNodes[1].push_back(numNodes + nb); polyTotNumNodes[1] += numNodes; if (numComp == 1){ NbSD++; return &SD; } else if(numComp == 3){ NbVD++; return &VD; } else if(numComp == 9){ NbTD++; return &TD; } break; } return 0; }