Write mesh to output file. 
   70{
   72 
   73    
   74    if (
m_f->m_exp[0]->GetNumElmts() == 0)
 
   75    {
   76        return;
   77    }
   78 
   80             "Need to specify fromfld=file.fld ");
   81 
   82    int nstart, i, j, nfields = 0;
   83    bool wssg      = false;
   85    string fromfld, basename, endname, nstartStr;
   86    stringstream filename;
   87    vector<string> infiles(nfld);
   88    vector<std::shared_ptr<Field>> fromField(nfld);
   89 
   90    
   91    fromfld  = 
m_config[
"fromfld"].as<
string>();
 
   92    basename = fromfld.substr(0, fromfld.find_first_of("_") + 1);
   93    filename << fromfld.substr(fromfld.find_first_of("_") + 1, fromfld.size());
   94    filename >> nstart;
   95    filename.str("");
   96    filename << nstart;
   97    filename >> nstartStr;
   98    filename.str("");
   99    endname = fromfld.substr(fromfld.find(nstartStr) + nstartStr.size(),
  100                             fromfld.size());
  101 
  102    for (i = 0; i < nfld; ++i)
  103    {
  104        stringstream filename;
  105        filename << basename << i + nstart << endname;
  106        filename >> infiles[i];
  107        cout << infiles[i] << endl;
  108    }
  109 
  110    for (i = 0; i < nfld; ++i)
  111    {
  112        fromField[i]            = std::shared_ptr<Field>(new Field());
  113        fromField[i]->m_session = 
m_f->m_session;
 
  114        fromField[i]->m_graph   = 
m_f->m_graph;
 
  115    }
  116 
  117    
  118    for (i = 0; i < nfld; ++i)
  119    {
  120        if (
m_f->m_exp.size())
 
  121        {
  122            
  123            Array<OneD, int> ElementGIDs(
m_f->m_exp[0]->GetExpSize());
 
  124            for (j = 0; j < 
m_f->m_exp[0]->GetExpSize(); ++j)
 
  125            {
  126                ElementGIDs[j] =
  127                    m_f->m_exp[0]->GetExp(j)->GetGeom()->GetGlobalID();
 
  128            }
  129            m_f->FieldIOForFile(infiles[i])
 
  130                ->Import(infiles[i], fromField[i]->m_fielddef,
  131                         fromField[i]->m_data,
  133        }
  134        else
  135        {
  136            m_f->FieldIOForFile(infiles[i])
 
  137                ->Import(infiles[i], fromField[i]->m_fielddef,
  138                         fromField[i]->m_data,
  140        }
  141 
  142        nfields = fromField[i]->m_fielddef[0]->m_fields.size();
  143        int NumHomogeneousDir =
  144            fromField[i]->m_fielddef[0]->m_numHomogeneousDir;
  145 
  146        if (nfields == 5 || nfields == 7)
  147        {
  148            wssg = true;
  149        }
  150 
  151        
  152        fromField[i]->m_graph->SetExpansionInfo(fromField[i]->m_fielddef);
  153 
  154        
  155        fromField[i]->m_exp.resize(nfields);
  156        fromField[i]->m_exp[0] =
  157            fromField[i]->SetUpFirstExpList(NumHomogeneousDir, true);
  158 
  159        for (j = 1; j < nfields; ++j)
  160        {
  161            fromField[i]->m_exp[j] = 
m_f->AppendExpList(NumHomogeneousDir);
 
  162        }
  163 
  164        for (j = 0; j < nfields; ++j)
  165        {
  166            for (int k = 0; k < fromField[i]->m_data.size(); ++k)
  167            {
  168                fromField[i]->m_exp[j]->ExtractDataToCoeffs(
  169                    fromField[i]->m_fielddef[k], fromField[i]->m_data[k],
  170                    fromField[i]->m_fielddef[k]->m_fields[j],
  171                    fromField[i]->m_exp[j]->UpdateCoeffs());
  172            }
  173            fromField[i]->m_exp[j]->BwdTrans(
  174                fromField[i]->m_exp[j]->GetCoeffs(),
  175                fromField[i]->m_exp[j]->UpdatePhys());
  176        }
  177    }
  178 
  179    int spacedim = 
m_f->m_graph->GetSpaceDimension();
 
  180    if ((fromField[0]->m_fielddef[0]->m_numHomogeneousDir) == 1 ||
  181        (fromField[0]->m_fielddef[0]->m_numHomogeneousDir) == 2)
  182    {
  183        spacedim = 3;
  184    }
  185 
  186    int nout = 5; 
  187    if (wssg)
  188    {
  189        nout = 6;
  190    }
  191 
  192    int npoints = fromField[0]->m_exp[0]->GetNpoints();
  193    Array<OneD, Array<OneD, NekDouble>> normTemporalMeanVec(spacedim),
  194        normCrossDir(spacedim), outfield(nout), dtemp(spacedim);
  195    Array<OneD, NekDouble> TemporalMeanMag(npoints, 0.0),
  196        DotProduct(npoints, 0.0), temp(npoints, 0.0);
  197 
  198    for (i = 0; i < spacedim; ++i)
  199    {
  200        normTemporalMeanVec[i] = Array<OneD, NekDouble>(npoints);
  201        normCrossDir[i]        = Array<OneD, NekDouble>(npoints);
  202        dtemp[i]               = Array<OneD, NekDouble>(npoints);
  205    }
  206 
  207    for (i = 0; i < nout; ++i)
  208    {
  209        outfield[i] = Array<OneD, NekDouble>(npoints, 0.0);
  210    }
  211 
  212    
  213    
  214    
  215    for (i = 0; i < nfld; ++i)
  216    {
  217        for (j = 0; j < spacedim; ++j)
  218        {
  219            Vmath::Vadd(npoints, fromField[i]->m_exp[j]->GetPhys(), 1,
 
  220                        normTemporalMeanVec[j], 1, normTemporalMeanVec[j], 1);
  221        }
  222    }
  223 
  224    for (i = 0; i < spacedim; ++i)
  225    {
  226        Vmath::Smul(npoints, 1.0 / nfld, normTemporalMeanVec[i], 1,
 
  227                    normTemporalMeanVec[i], 1);
  228        Vmath::Vvtvp(npoints, normTemporalMeanVec[i], 1, normTemporalMeanVec[i],
 
  229                     1, TemporalMeanMag, 1, TemporalMeanMag, 1);
  230    }
  231 
  232    Vmath::Vsqrt(npoints, TemporalMeanMag, 1, TemporalMeanMag, 1);
 
  233 
  234    for (i = 0; i < spacedim; ++i)
  235    {
  236        Vmath::Vdiv(npoints, normTemporalMeanVec[i], 1, TemporalMeanMag, 1,
 
  237                    normTemporalMeanVec[i], 1);
  238    }
  239    
  240 
  241    if (wssg) 
  242              
  243    {
  244        Vmath::Vmul(npoints, fromField[0]->m_exp[nfields - 1]->GetPhys(), 1,
 
  245                    normTemporalMeanVec[1], 1, normCrossDir[0], 1);
  246        Vmath::Vvtvm(npoints, fromField[0]->m_exp[nfields - 2]->GetPhys(), 1,
 
  247                     normTemporalMeanVec[2], 1, normCrossDir[0], 1,
  248                     normCrossDir[0], 1);
  249        Vmath::Vmul(npoints, fromField[0]->m_exp[nfields - 3]->GetPhys(), 1,
 
  250                    normTemporalMeanVec[2], 1, normCrossDir[1], 1);
  251        Vmath::Vvtvm(npoints, fromField[0]->m_exp[nfields - 1]->GetPhys(), 1,
 
  252                     normTemporalMeanVec[0], 1, normCrossDir[1], 1,
  253                     normCrossDir[1], 1);
  254        Vmath::Vmul(npoints, fromField[0]->m_exp[nfields - 2]->GetPhys(), 1,
 
  255                    normTemporalMeanVec[0], 1, normCrossDir[2], 1);
  256        Vmath::Vvtvm(npoints, fromField[0]->m_exp[nfields - 3]->GetPhys(), 1,
 
  257                     normTemporalMeanVec[1], 1, normCrossDir[2], 1,
  258                     normCrossDir[2], 1);
  259    }
  260 
  261    
  262    for (i = 0; i < nfld; ++i)
  263    {
  264        for (j = 0; j < spacedim; ++j)
  265        {
  266            Vmath::Vvtvp(npoints, fromField[i]->m_exp[j]->GetPhys(), 1,
 
  267                         normTemporalMeanVec[j], 1, DotProduct, 1, DotProduct,
  268                         1);
  269        }
  270        
  271        Vmath::Vadd(npoints, fromField[i]->m_exp[spacedim]->GetPhys(), 1,
 
  272                    outfield[0], 1, outfield[0], 1);
  273        
  274        Vmath::Vmul(npoints, DotProduct, 1, DotProduct, 1, temp, 1);
 
  275        Vmath::Vvtvm(npoints, fromField[i]->m_exp[spacedim]->GetPhys(), 1,
 
  276                     fromField[i]->m_exp[spacedim]->GetPhys(), 1, temp, 1, temp,
  277                     1);
  278 
  279        for (j = 0; j < npoints; ++j)
  280        {
  281            if (temp[j] > 0.0)
  282            {
  283                outfield[1][j] = outfield[1][j] + 
sqrt(temp[j]);
 
  284            }
  285        }
  286 
  287        
  289                    fromField[i]->m_exp[spacedim]->GetPhys(), 1, temp, 1);
  290        Vmath::Vadd(npoints, temp, 1, outfield[3], 1, outfield[3], 1);
 
  291 
  292        
  293        for (j = 0; j < npoints; ++j)
  294        {
  295            temp[j] = 1 - temp[j] * temp[j];
  296            if (temp[j] > 0.0)
  297            {
  298                outfield[4][j] = outfield[4][j] + 
sqrt(temp[j]);
 
  299            }
  300        }
  301 
  302        
  303        if (wssg)
  304        {
  305            
  307 
  308            fromField[i]->m_exp[0]->PhysDeriv(DotProduct, dtemp[0], dtemp[1],
  309                                              dtemp[2]);
  310            for (j = 0; j < spacedim; j++)
  311            {
  312                Vmath::Vvtvp(npoints, dtemp[j], 1, normTemporalMeanVec[j], 1,
 
  313                             temp, 1, temp, 1);
  314            }
  316 
  317            
  319 
  320            for (j = 0; j < spacedim; ++j)
  321            {
  322                Vmath::Vvtvp(npoints, fromField[i]->m_exp[j]->GetPhys(), 1,
 
  323                             normCrossDir[j], 1, DotProduct, 1, DotProduct, 1);
  324            }
  325            fromField[i]->m_exp[0]->PhysDeriv(DotProduct, dtemp[0], dtemp[1],
  326                                              dtemp[2]);
  328 
  329            for (j = 0; j < spacedim; j++)
  330            {
  332                             DotProduct, 1, DotProduct, 1);
  333            }
  334            Vmath::Vvtvp(npoints, DotProduct, 1, DotProduct, 1, temp, 1, temp,
 
  335                         1);
  336 
  337            
  339            Vmath::Vadd(npoints, temp, 1, outfield[5], 1, outfield[5], 1);
 
  340        }
  341 
  343    }
  344 
  345    
  346    Vmath::Smul(npoints, 1.0 / nfld, outfield[0], 1, outfield[0], 1);
 
  347    Vmath::Smul(npoints, 1.0 / nfld, outfield[1], 1, outfield[1], 1);
 
  348    Vmath::Smul(npoints, 1.0 / nfld, outfield[3], 1, outfield[3], 1);
 
  349    Vmath::Smul(npoints, 1.0 / nfld, outfield[4], 1, outfield[4], 1);
 
  350 
  351    
  352    for (i = 0; i < npoints; ++i)
  353    {
  354        outfield[2][i] = 0.5 * (1 - TemporalMeanMag[i] / outfield[0][i]);
  355    }
  356 
  357    
  358
  359
  360
  361
  362
  363 
  364    m_f->m_exp.resize(nout);
 
  365    m_f->m_fielddef = fromField[0]->m_fielddef;
 
  367        m_f->SetUpFirstExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir, 
true);
 
  368 
  369    for (i = 1; i < nout; ++i)
  370    {
  372            m_f->AppendExpList(
m_f->m_fielddef[0]->m_numHomogeneousDir);
 
  373    }
  374 
  375    m_f->m_fielddef[0]->m_fields.resize(nout);
 
  376    m_f->m_fielddef[0]->m_fields[0] = 
"TAWSS";
 
  377    m_f->m_fielddef[0]->m_fields[1] = 
"transWSS";
 
  378    m_f->m_fielddef[0]->m_fields[2] = 
"OSI";
 
  379    m_f->m_fielddef[0]->m_fields[3] = 
"TAAFI";
 
  380    m_f->m_fielddef[0]->m_fields[4] = 
"TACFI";
 
  381 
  382    if (wssg)
  383    {
  384        m_f->m_fielddef[0]->m_fields[5] = 
"|WSSG|";
 
  385        Vmath::Smul(npoints, 1.0 / nfld, outfield[5], 1, outfield[5], 1);
 
  386    }
  387 
  388    m_f->m_variables = 
m_f->m_fielddef[0]->m_fields;
 
  389 
  390    for (i = 0; i < nout; ++i)
  391    {
  392        m_f->m_exp[i]->FwdTrans(outfield[i], 
m_f->m_exp[i]->UpdateCoeffs());
 
  393        m_f->m_exp[i]->BwdTrans(
m_f->m_exp[i]->GetCoeffs(),
 
  394                                m_f->m_exp[i]->UpdatePhys());
 
  395    }
  396}
#define ASSERTL0(condition, msg)
FieldSharedPtr m_f
Field object.
static FieldMetaDataMap NullFieldMetaDataMap
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Vvtvm(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvm (vector times vector minus vector): z = w*x - y
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
void Zero(int n, T *x, const int incx)
Zero vector.
scalarT< T > sqrt(scalarT< T > in)