← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to lib/atc/FE_Engine.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/atc/FE_Engine.cpp
89
89
    // now do all FE_Engine data structure partitioning
90
90
 
91
91
    // partition nullElements_
92
 
    /*for (vector<int>::const_iterator elemsIter = feMesh_->processor_elts().begin();
 
92
    /*for (vector<int>::iterator elemsIter = feMesh_->processor_elts().begin();
93
93
         elemsIter != feMesh_->processor_elts().end();
94
94
         ++elemsIter)
95
95
    {
182
182
          dx = 0;
183
183
          dy = 0;
184
184
          dz = 0;
185
 
          double x[3] = {0,0,0}; 
186
185
          while (argIdx < narg) {
187
 
            if      (strcmp(arg[argIdx++],"dx")==0) {
188
 
              // parse relative values for each element
189
 
              if (is_numeric(arg[argIdx])) {
190
 
                for (int i = 0; i < nx; ++i) { 
191
 
                  if (is_numeric(arg[argIdx])) { dx(i) = atof(arg[argIdx++]); }
192
 
                  else throw ATC_Error("not enough element partitions");
193
 
                }
194
 
              }
195
 
              // construct relative values from a density function
196
 
              // evaluate for a domain (0,1)
197
 
              else {
198
 
                XT_Function * f = XT_Function_Mgr::instance()->function(&(arg[argIdx]),narg-argIdx);
199
 
                for (int i = 0; i < nx; ++i) { 
200
 
                  x[0] = (i+0.5)/nx; dx(i) = f->f(x,0.); 
201
 
                }
202
 
              }
203
 
 
204
 
 
205
 
              break ;
206
 
            }
207
 
            else if (strcmp(arg[argIdx++],"dy")==0) {
208
 
              for (int i = 0; i < ny; ++i) { dy(i) = atof(arg[argIdx++]); }
209
 
            }
210
 
            else if (strcmp(arg[argIdx++],"dz")==0) {
211
 
              for (int i = 0; i < nz; ++i) { dz(i) = atof(arg[argIdx++]); }
 
186
            if      (strcmp(arg[argIdx],"dx")==0) {
 
187
              parse_partitions(++argIdx, narg, arg, 0, dx);
 
188
            }
 
189
            else if (strcmp(arg[argIdx],"dy")==0) {
 
190
              parse_partitions(++argIdx, narg, arg, 1, dy);
 
191
            }
 
192
            else if (strcmp(arg[argIdx],"dz")==0) {
 
193
              parse_partitions(++argIdx, narg, arg, 2, dz);
212
194
            }
213
195
          }
214
196
 
354
336
    }
355
337
    return match;
356
338
  }
 
339
  //-----------------------------------------------------------------
 
340
  void FE_Engine::parse_partitions(int & argIdx, int narg, char ** arg, 
 
341
    int idof, Array<double> & dx ) const
 
342
  {
 
343
    double x[3] = {0,0,0}; 
 
344
    int nx = dx.size();
 
345
    // parse relative values for each element
 
346
    if (is_numeric(arg[argIdx])) {
 
347
      for (int i = 0; i < nx; ++i) { 
 
348
        if (is_numeric(arg[argIdx])) { dx(i) = atof(arg[argIdx++]); }
 
349
        else throw ATC_Error("not enough element partitions");
 
350
      }
 
351
    }
 
352
    // each segment of the piecewise funcion is length-normalized separately
 
353
    else if (strcmp(arg[argIdx],"position-number-density")==0) { 
 
354
      argIdx++;
 
355
      double y[nx],w[nx];
 
356
      int n[nx];
 
357
      int nn = 0;
 
358
      while (argIdx < narg) { 
 
359
        if (! is_numeric(arg[argIdx])) break;
 
360
        y[nn]   = atof(arg[argIdx++]);
 
361
        n[nn]   = atoi(arg[argIdx++]);
 
362
        w[nn++] = atof(arg[argIdx++]);
 
363
      }
 
364
      if (n[nn-1] != nx)  throw ATC_Error("total element partitions do not match");
 
365
      int k = 0;
 
366
      for (int i = 1; i < nn; ++i) { 
 
367
        int dn = n[i]-n[i-1];
 
368
        double dy = y[i]-y[i-1];
 
369
        double w0 = w[i-1];
 
370
        double dw = w[i]-w0;
 
371
        double lx = 0;
 
372
        double l[dn];
 
373
        for (int j = 0; j < dn; ++j) {
 
374
          double x = (j+0.5)/dn; 
 
375
          double dl = w0+x*dw;
 
376
          lx += dl;
 
377
          l[j]= dl;
 
378
        }
 
379
        double scale = dy/lx;
 
380
        for (int j = 0; j < dn; ++j) {
 
381
          dx(k++) = scale*l[j];
 
382
        }
 
383
      }
 
384
    }
 
385
    // construct relative values from a density function
 
386
    // evaluate for a domain (0,1)
 
387
    else {
 
388
      XT_Function * f = XT_Function_Mgr::instance()->function(&(arg[argIdx]),narg-argIdx);
 
389
      argIdx++;
 
390
      while (argIdx < narg) { 
 
391
        if (! is_numeric(arg[argIdx++])) break;
 
392
      }
 
393
      for (int i = 0; i < nx; ++i) { 
 
394
        x[idof] = (i+0.5)/nx; dx(i) = f->f(x,0.); 
 
395
      }
 
396
    }
 
397
  }
357
398
 
358
399
  //-----------------------------------------------------------------
359
400
  void FE_Engine::finish()
550
591
    DENS_MAT elementMassMatrix(nNodesPerElement_,nNodesPerElement_); 
551
592
 
552
593
    vector<int> myElems = feMesh_->owned_elts();
553
 
    for (vector<int>::const_iterator elemsIter = myElems.begin();
 
594
    for (vector<int>::iterator elemsIter = myElems.begin();
554
595
         elemsIter != myElems.end();
555
596
         ++elemsIter)
556
597
    {
604
645
    DENS_MAT coefsAtIPs;
605
646
 
606
647
    vector<int> myElems = feMesh_->owned_elts();
607
 
    for (vector<int>::const_iterator elemsIter = myElems.begin();
 
648
    for (vector<int>::iterator elemsIter = myElems.begin();
608
649
         elemsIter != myElems.end();
609
650
         ++elemsIter)
610
651
    {
782
823
 
783
824
    // element wise assembly
784
825
    vector<int> myElems = feMesh_->owned_elts();
785
 
    for (vector<int>::const_iterator elemsIter = myElems.begin(); 
 
826
    for (vector<int>::iterator elemsIter = myElems.begin(); 
786
827
         elemsIter != myElems.end();
787
828
         ++elemsIter) 
788
829
    {
849
890
    massMatrix.reset(nNodesUnique_,nNodesUnique_);// zero since partial fill
850
891
    DENS_MAT elementMassMatrix(nNodesPerElement_,nNodesPerElement_); 
851
892
    vector<int> myElems = feMesh_->owned_elts();
852
 
    for (vector<int>::const_iterator elemsIter = myElems.begin();
 
893
    for (vector<int>::iterator elemsIter = myElems.begin();
853
894
         elemsIter != myElems.end();
854
895
         ++elemsIter)
855
896
    {
905
946
    // assemble lumped diagonal mass 
906
947
    DENS_VEC Nvec(nNodesPerElement_);
907
948
    vector<int> myElems = feMesh_->owned_elts();
908
 
    for (vector<int>::const_iterator elemsIter = myElems.begin();
 
949
    for (vector<int>::iterator elemsIter = myElems.begin();
909
950
         elemsIter != myElems.end();
910
951
         ++elemsIter)
911
952
    {
944
985
    }
945
986
    // assemble diagonal matrix
946
987
    vector<int> myElems = feMesh_->owned_elts();
947
 
    for (vector<int>::const_iterator elemsIter = myElems.begin();
 
988
    for (vector<int>::iterator elemsIter = myElems.begin();
948
989
         elemsIter != myElems.end();
949
990
         ++elemsIter)
950
991
    {
1028
1069
      // treat single material point sets specially
1029
1070
      int nMatls = pointMaterialGroups.size();
1030
1071
      int atomMatls = 0;
 
1072
      int iatomMat = 0;
1031
1073
      for (int imat = 0; imat < nMatls; imat++) {
1032
1074
        const set<int> & indices = pointMaterialGroups(imat);
1033
 
        if ( indices.size() > 0) atomMatls++;
 
1075
        if ( indices.size() > 0) {
 
1076
          iatomMat = imat;
 
1077
          atomMatls++;
 
1078
        }
1034
1079
      }
1035
1080
      if (atomMatls == 0) {
1036
1081
        throw ATC_Error("no materials in atom region - atoms may have migrated to FE-only region");
1045
1090
      FIELD_MATS capacities;
1046
1091
      // evaluate physics model & compute capacity|density for requested fields
1047
1092
      if (singleMaterial) {
1048
 
        int imat = 0; 
1049
 
        const Material * mat = physicsModel->material(imat);
 
1093
        const Material * mat = physicsModel->material(iatomMat);
1050
1094
        for (int j = 0; j < nfields; ++j) {
1051
1095
          _fieldName_ = field_mask(j);
1052
1096
          // mass matrix needs to be invertible so null matls have cap=1
1053
 
          if (physicsModel->null(_fieldName_,imat)) {
 
1097
          if (physicsModel->null(_fieldName_,iatomMat)) {
1054
1098
             throw ATC_Error("null material not supported for atomic region (single material)");
1055
1099
             const FIELD &  f = (fields.find(_fieldName_))->second;
1056
1100
             capacities[_fieldName_].reset(f.nRows(),f.nCols());
1143
1187
    Array<int> count(nNodesUnique_); count = 0;
1144
1188
    set_quadrature(NODAL);
1145
1189
    vector<int> myElems = feMesh_->owned_elts();
1146
 
    for (vector<int>::const_iterator elemsIter = myElems.begin();
 
1190
    for (vector<int>::iterator elemsIter = myElems.begin();
1147
1191
         elemsIter != myElems.end();
1148
1192
         ++elemsIter)
1149
1193
    {
1204
1248
 
1205
1249
    DENS_MAT elementEnergy(nNodesPerElement_,1); // workspace
1206
1250
    vector<int> myElems = feMesh_->owned_elts();
1207
 
    for (vector<int>::const_iterator elemsIter = myElems.begin();
 
1251
    for (vector<int>::iterator elemsIter = myElems.begin();
1208
1252
         elemsIter != myElems.end();
1209
1253
         ++elemsIter)
1210
1254
    {
1249
1293
    const PhysicsModel * physicsModel,
1250
1294
    const Array<int>   & elementMaterials,
1251
1295
    FIELDS &rhs,
 
1296
    bool freeOnly,
1252
1297
    const DenseMatrix<bool> *elementMask) const
1253
1298
  {
1254
1299
    vector<FieldName> usedFields;
1269
1314
 
1270
1315
    // Iterate over elements partitioned onto current processor.
1271
1316
    vector<int> myElems = feMesh_->owned_elts();
1272
 
    for (vector<int>::const_iterator elemsIter = myElems.begin(); 
 
1317
    for (vector<int>::iterator elemsIter = myElems.begin(); 
1273
1318
         elemsIter != myElems.end();
1274
1319
         ++elemsIter) 
1275
1320
    {
1337
1382
      } 
1338
1383
    }
1339
1384
 
1340
 
    vector<FieldName>::const_iterator _fieldIter_;
 
1385
    vector<FieldName>::iterator _fieldIter_;
1341
1386
    for (_fieldIter_ = usedFields.begin(); _fieldIter_ != usedFields.end(); 
1342
1387
         ++_fieldIter_) {
1343
1388
      // myRhs is where we put the global result for this field.
1660
1705
 
1661
1706
    // element wise assembly
1662
1707
    vector<int> myElems = feMesh_->owned_elts();
1663
 
    for (vector<int>::const_iterator elemsIter = myElems.begin();
 
1708
    for (vector<int>::iterator elemsIter = myElems.begin();
1664
1709
         elemsIter != myElems.end();
1665
1710
         ++elemsIter)
1666
1711
    {
1740
1785
    }
1741
1786
  
1742
1787
    // element wise assembly
1743
 
    set< PAIR >::const_iterator iter;
 
1788
    set< PAIR >::iterator iter;
1744
1789
    for (iter = faceSet.begin(); iter != faceSet.end(); iter++) {
1745
1790
      // get connectivity
1746
1791
      int ielem = iter->first;
1841
1886
    
1842
1887
    // R_I    = - int_Omega    Delta _N_I .q dV
1843
1888
    compute_rhs_vector(rhsMask, fields, physicsModel, elementMaterials, rhs, elementMask);
 
1889
 
1844
1890
    // R_I^md = - int_Omega^md Delta _N_I .q dV
1845
1891
    compute_rhs_vector(rhsMask, fields, physicsModel, pointMaterialGroups, 
1846
1892
      _weights_, N, dN, rhsSubset);
1847
 
    
 
1893
 
1848
1894
    // flux_I = 1/Delta V_I V_I^md R_I + R_I^md
1849
1895
    //   where : Delta V_I = int_Omega _N_I dV
1850
1896
    for (int n = 0; n < rhsMask.nRows(); n++) {
1920
1966
    DENS_MAT localField;
1921
1967
    DENS_MAT xAtIPs(nSD_,nIPsPerFace_);
1922
1968
    DENS_MAT uAtIPs(nIPsPerFace_,1); 
1923
 
 
 
1969
    FIELDS myNodalSources;
1924
1970
 
1925
1971
    // element wise assembly
1926
1972
    ROBIN_SURFACE_SOURCE::const_iterator src_iter;
1927
 
    if (!(rank_zero(communicator_))) {
1928
 
      // Zero out unmasked nodal sources on all non-main processors.
1929
 
      // This is to avoid counting the previous nodal source values
1930
 
      // multiple times when we aggregate.
1931
 
      for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
1932
 
      {
1933
 
        _fieldName_ = src_iter->first;
1934
 
        if (!rhsMask((int)_fieldName_,ROBIN_SOURCE)) continue; 
1935
 
        DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
1936
 
        myNodalSources.reset(myNodalSources.nRows(), myNodalSources.nCols());
1937
 
      }
 
1973
    for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
 
1974
    {
 
1975
      _fieldName_ = src_iter->first;
 
1976
      if (!rhsMask((int)_fieldName_,ROBIN_SOURCE)) continue; 
 
1977
      DENS_MAT & s(nodalSources[_fieldName_].set_quantity());
 
1978
      myNodalSources[_fieldName_] = DENS_MAT(s.nRows(),s.nCols());
1938
1979
    }
1939
1980
    for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
1940
1981
    {
1982
2023
        }
1983
2024
        
1984
2025
        // assemble
1985
 
        DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
 
2026
        DENS_MAT & s(myNodalSources[_fieldName_].set_quantity());
1986
2027
        _Nmat_ = _fN_.transMat(_fweights_*faceSource);
1987
2028
        for (int i = 0; i < nNodesPerElement_; ++i) {
1988
2029
          int inode = _conn_(i);
1989
2030
          for (int idof = 0; idof < nFieldDOF; ++idof) {
1990
 
            myNodalSources(inode,idof) += _Nmat_(i,idof);
 
2031
            s(inode,idof) += _Nmat_(i,idof);
1991
2032
          }  
1992
2033
        } 
1993
2034
      }  
1996
2037
    for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++) {
1997
2038
      _fieldName_ = src_iter->first;
1998
2039
      if (!rhsMask((int) _fieldName_,ROBIN_SOURCE)) continue; 
1999
 
      DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
2000
 
      allsum(communicator_,MPI_IN_PLACE, myNodalSources.ptr(), myNodalSources.size());  
 
2040
      DENS_MAT & s(myNodalSources[_fieldName_].set_quantity());
 
2041
      allsum(communicator_,MPI_IN_PLACE, s.ptr(), s.size());  
 
2042
      DENS_MAT & src(nodalSources[_fieldName_].set_quantity());
 
2043
      src += s;
2001
2044
    }
2002
2045
  }
2003
2046
 
2021
2064
 
2022
2065
    // element wise assembly
2023
2066
    ROBIN_SURFACE_SOURCE::const_iterator src_iter;
2024
 
    if (!(rank_zero(communicator_))) {
2025
 
      // Zero out result (tangent) matrix on all non-main processors
2026
 
      // to avoid multiple-counting of the values already in tangent 
2027
 
      tangent.reset(tangent.nRows(), tangent.nCols());
2028
 
    }
 
2067
    SPAR_MAT K(tangent.nRows(), tangent.nCols());
2029
2068
    for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
2030
2069
    {
2031
2070
      _fieldName_ = src_iter->first;
2068
2107
        // assemble
2069
2108
        DIAG_MAT D(column(coefsAtIPs,0)); 
2070
2109
        D *= -1; 
2071
 
 
2072
2110
        D *= _fweights_; 
2073
2111
        _Nmat_ = _fN_.transMat(D*_fN_);
2074
2112
        for (int i = 0; i < nNodesPerElement_; ++i) {
2075
2113
          int inode = _conn_(i);
2076
2114
          for (int j = 0; j < nNodesPerElement_; ++j) {
2077
2115
            int jnode = _conn_(j);
2078
 
            tangent.add(inode, jnode, _Nmat_(i,j));
2079
 
          }  
2080
 
        } 
2081
 
      }  
2082
 
    } 
2083
 
    // assemble partial result matrices
2084
 
#ifdef ISOLATE_FE
2085
 
    sparse_allsum(communicator_,tangent);
2086
 
#else
2087
 
    LammpsInterface::instance()->sparse_allsum(tangent);
2088
 
#endif
 
2116
            K.add(inode, jnode, _Nmat_(i,j));
 
2117
          }  
 
2118
        } 
 
2119
      }  
 
2120
    } 
 
2121
    // assemble partial result matrices
 
2122
#ifdef ISOLATE_FE
 
2123
    sparse_allsum(communicator_,K);
 
2124
#else
 
2125
    LammpsInterface::instance()->sparse_allsum(K);
 
2126
#endif
 
2127
    tangent += K;
 
2128
  }
 
2129
 
 
2130
  //-----------------------------------------------------------------
 
2131
  // Robin boundary flux using native quadrature 
 
2132
  // integrate \int_delV _N_I s(x,t).n dA
 
2133
  //-----------------------------------------------------------------
 
2134
  void FE_Engine::add_open_fluxes(const Array2D<bool> &rhsMask,
 
2135
                                  const FIELDS        & fields,
 
2136
                                  const OPEN_SURFACE  & openFaces,
 
2137
                                  FIELDS &nodalSources,
 
2138
                                  const FieldName Velocity) const
 
2139
  {
 
2140
    
 
2141
    // sizing working arrays
 
2142
    DENS_MAT xCoords(nSD_,nNodesPerElement_);
 
2143
    DENS_MAT faceSource;
 
2144
    DENS_MAT localField;
 
2145
    DENS_MAT xAtIPs(nSD_,nIPsPerFace_);
 
2146
    DENS_MAT uAtIPs; // (nIPsPerFace_,field nCols);
 
2147
    DENS_MAT aAtIPs(nIPsPerFace_,1);
 
2148
    FIELDS myNodalSources;
 
2149
 
 
2150
    // throw error if electron velocity is not defined
 
2151
 
 
2152
    // element wise assembly
 
2153
    OPEN_SURFACE::const_iterator face_iter;
 
2154
    for (face_iter=openFaces.begin(); face_iter!=openFaces.end(); face_iter++)
 
2155
    {
 
2156
      _fieldName_ = face_iter->first;
 
2157
      if (!rhsMask((int)_fieldName_,OPEN_SOURCE)) continue; 
 
2158
      DENS_MAT & s(nodalSources[_fieldName_].set_quantity());
 
2159
      myNodalSources[_fieldName_] = DENS_MAT(s.nRows(),s.nCols());
 
2160
    }
 
2161
    for (face_iter=openFaces.begin(); face_iter!=openFaces.end(); face_iter++)
 
2162
    {
 
2163
      _fieldName_ = face_iter->first;
 
2164
      if (!rhsMask((int)_fieldName_,OPEN_SOURCE)) continue; 
 
2165
      
 
2166
      typedef set<PAIR> FSET;
 
2167
      const FSET *fset = (const FSET *)&(face_iter->second);
 
2168
      FSET::const_iterator fset_iter;
 
2169
      for (fset_iter = fset->begin(); fset_iter != fset->end(); fset_iter++)
 
2170
      {
 
2171
        const PAIR &face = *fset_iter;
 
2172
        const int elem = face.first;
 
2173
        // if this is not our element, do not do calculations
 
2174
        if (!feMesh_->is_owned_elt(elem)) continue;
 
2175
        // get velocity, multiply with field (vector gives tensor)
 
2176
        // dot with face normal
 
2177
        _conn_ = feMesh_->element_connectivity_unique(elem);
 
2178
        // evaluate location at ips
 
2179
        feMesh_->face_shape_function(face, _fN_, _fdN_, _nN_, _fweights_);
 
2180
        // collect field at ips
 
2181
        _fieldItr_ = fields.find(_fieldName_);
 
2182
        const DENS_MAT & field = (_fieldItr_->second).quantity();
 
2183
        feMesh_->element_field(elem, field, localField);
 
2184
        int nFieldDOF = field.nCols();
 
2185
        // "u" is the quantity being advected
 
2186
        uAtIPs.reset(nIPsPerFace_,nFieldDOF);
 
2187
        uAtIPs = _fN_*localField;
 
2188
        // collect velocity at ips for "advection" = v.n
 
2189
        _fieldItr_ = fields.find(Velocity);
 
2190
        const DENS_MAT & advection = (_fieldItr_->second).quantity(); // advection velocity
 
2191
        feMesh_->element_field(elem, advection, localField);
 
2192
        for (int j = 0; j < nSD_; ++j) { // nSD_ == nDOF for the velocity
 
2193
          for (int ip = 0; ip < nIPsPerFace_; ++ip) {
 
2194
            for (int I = 0; I < nNodesPerElement_; ++I) {
 
2195
              aAtIPs(ip,0) += (_nN_[j])(ip,I)*localField(I,j);
 
2196
            }
 
2197
          }
 
2198
        }
 
2199
        // construct open flux at ips of this element
 
2200
        // flux = field \otimes advection_vector \dot n
 
2201
        faceSource.reset(nIPsPerFace_,nFieldDOF);
 
2202
        for (int ip = 0; ip < nIPsPerFace_; ++ip) {
 
2203
          for (int idof = 0; idof<nFieldDOF; ++idof) {
 
2204
            // multiply field DOF value times advection vector
 
2205
            faceSource(ip,idof) = aAtIPs(ip,1)*uAtIPs(ip,idof);//(v.n) u
 
2206
          }
 
2207
        }
 
2208
        // assemble contributions for this direction of the face normal
 
2209
        // _nN_[j](ip,I) nodal shapefunction(I) at ip X face normal(j)
 
2210
        // _fweights_ diagonal nips X nips ==> facesource is nips X ndofs
 
2211
        DENS_MAT & s(myNodalSources[_fieldName_].set_quantity());
 
2212
        _Nmat_ = _fN_.transMat(_fweights_*faceSource);
 
2213
        // s_Ii = \sum_ip N_I (u_i v.n)_ip wg_ip
 
2214
        for (int i = 0; i < nNodesPerElement_; ++i) {
 
2215
          int inode = _conn_(i);
 
2216
          for (int idof = 0; idof < nFieldDOF; ++idof) {
 
2217
            s(inode,idof) += _Nmat_(i,idof);
 
2218
          }  
 
2219
        }
 
2220
      }  
 
2221
    } 
 
2222
    // assemble partial result matrices
 
2223
    for (face_iter=openFaces.begin(); face_iter!=openFaces.end(); face_iter++) {
 
2224
      _fieldName_ = face_iter->first;
 
2225
      if (!rhsMask((int) _fieldName_,OPEN_SOURCE)) continue; 
 
2226
      DENS_MAT & s(myNodalSources[_fieldName_].set_quantity());
 
2227
      allsum(communicator_,MPI_IN_PLACE, s.ptr(), s.size());  
 
2228
      DENS_MAT & src(nodalSources[_fieldName_].set_quantity());
 
2229
      src += s;
 
2230
    }
 
2231
  }
 
2232
 
 
2233
  //-----------------------------------------------------------------
 
2234
  // Open boundary flux stiffness using native quadrature 
 
2235
  // integrate \int_delV _N_I ds/du(x,t).n dA
 
2236
  //-----------------------------------------------------------------
 
2237
  void FE_Engine::add_open_tangent(const Array2D<bool> &rhsMask,
 
2238
    const FIELDS        & fields,
 
2239
    const OPEN_SURFACE  & openFaces,
 
2240
    SPAR_MAT &tangent,
 
2241
    const FieldName Velocity) const
 
2242
  {
 
2243
    
 
2244
    // sizing working arrays
 
2245
    DENS_MAT xCoords(nSD_,nNodesPerElement_);
 
2246
    DENS_MAT faceSource;
 
2247
    DENS_MAT localField;
 
2248
    DENS_MAT xAtIPs(nSD_,nIPsPerFace_);
 
2249
    DENS_MAT uAtIPs; // (nIPsPerFace_,field nCols);
 
2250
    DENS_MAT aAtIPs(nIPsPerFace_,nSD_);
 
2251
    SPAR_MAT K(tangent.nRows(), tangent.nCols());
 
2252
    // element wise assembly
 
2253
    OPEN_SURFACE::const_iterator face_iter;
 
2254
    for (face_iter=openFaces.begin(); face_iter!=openFaces.end(); face_iter++)
 
2255
    {
 
2256
      _fieldName_ = face_iter->first;
 
2257
      if (!rhsMask((int)_fieldName_,OPEN_SOURCE)) continue; 
 
2258
      bool convective = false;
 
2259
      if (_fieldName_ == Velocity) convective = true;
 
2260
      
 
2261
      typedef set<PAIR> FSET;
 
2262
      const FSET *fset = (const FSET *)&(face_iter->second);
 
2263
      FSET::const_iterator fset_iter;
 
2264
      for (fset_iter = fset->begin(); fset_iter != fset->end(); fset_iter++)
 
2265
      {
 
2266
        const PAIR &face = *fset_iter;
 
2267
        const int elem = face.first;
 
2268
        // if this is not our element, do not do calculations
 
2269
        if (!feMesh_->is_owned_elt(elem)) continue;
 
2270
        _conn_ = feMesh_->element_connectivity_unique(elem);
 
2271
        // evaluate location at ips
 
2272
        feMesh_->face_shape_function(face, _fN_, _fdN_, _nN_, _fweights_);
 
2273
        // collect field at ips
 
2274
        _fieldItr_ = fields.find(_fieldName_);
 
2275
        const DENS_MAT & field = (_fieldItr_->second).quantity();
 
2276
        feMesh_->element_field(elem, field, localField);
 
2277
        int nFieldDOF = field.nCols();
 
2278
        // "u" is the quantity being advected
 
2279
        uAtIPs.reset(nIPsPerFace_,nFieldDOF);
 
2280
        uAtIPs = _fN_*localField;
 
2281
        // collect velocity at ips for "advection" = v.n
 
2282
        _fieldItr_ = fields.find(Velocity);
 
2283
        const DENS_MAT & advection = (_fieldItr_->second).quantity(); // advection velocity
 
2284
        feMesh_->element_field(elem, advection, localField);
 
2285
        for (int j = 0; j < nSD_; ++j) { // nSD_ == nDOF for the velocity
 
2286
          for (int ip = 0; ip < nIPsPerFace_; ++ip) {
 
2287
            for (int I = 0; I < nNodesPerElement_; ++I) {
 
2288
              aAtIPs(ip,0) += (_nN_[j])(ip,I)*localField(I,j);
 
2289
            }
 
2290
          }
 
2291
        }
 
2292
        // K_IJ = \sum_ip N_I ( v.n I + u (x) n ) N_J wg_ip
 
2293
        DENS_MAT D(nFieldDOF,nFieldDOF);
 
2294
        DENS_VEC n(nSD_);
 
2295
        for (int ip = 0; ip < nIPsPerFace_; ++ip) {
 
2296
          for (int idof = 0; idof<nFieldDOF; ++idof) {
 
2297
            D(idof,idof) -= aAtIPs(ip,0);
 
2298
            if (convective) {
 
2299
              feMesh_->face_normal(face,ip,n);
 
2300
              for (int jdof = 0; jdof<nFieldDOF; ++jdof) {
 
2301
                D(idof,jdof) -=  uAtIPs(ip,idof)*n(jdof);
 
2302
              }
 
2303
            }
 
2304
          }
 
2305
        }
 
2306
        // assemble
 
2307
        _Nmat_ = _fN_.transMat(D*_fN_);
 
2308
        for (int i = 0; i < nNodesPerElement_; ++i) {
 
2309
          int inode = _conn_(i);
 
2310
          for (int j = 0; j < nNodesPerElement_; ++j) {
 
2311
            int jnode = _conn_(j);
 
2312
            K.add(inode, jnode, _Nmat_(i,j));
 
2313
          }  
 
2314
        } 
 
2315
      }  
 
2316
    } 
 
2317
    // assemble partial result matrices
 
2318
#ifdef ISOLATE_FE
 
2319
    sparse_allsum(communicator_,K);
 
2320
#else
 
2321
    LammpsInterface::instance()->sparse_allsum(K);
 
2322
#endif
 
2323
    tangent += K;
2089
2324
  }
2090
2325
  
2091
2326
  //-----------------------------------------------------------------
2102
2337
    DENS_MAT xCoords(nSD_,nNodesPerElement_);
2103
2338
    DENS_MAT xAtIPs(nSD_,nIPsPerFace_);
2104
2339
    DENS_MAT faceSource;
 
2340
    FIELDS myNodalSources;
2105
2341
 
2106
2342
    // element wise assembly
2107
2343
    SURFACE_SOURCE::const_iterator src_iter;
2108
 
    if (!(rank_zero(communicator_))) {
2109
 
      // Zero out unmasked nodal sources on all non-main processors.
2110
 
      // This is to avoid counting the previous nodal source values
2111
 
      // multiple times when we aggregate.
2112
 
      for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
2113
 
      {
2114
 
        _fieldName_ = src_iter->first;
2115
 
        if (!fieldMask((int)_fieldName_)) continue; 
2116
 
        DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
2117
 
        myNodalSources.reset(myNodalSources.nRows(), myNodalSources.nCols());
2118
 
      }
 
2344
    for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
 
2345
    {
 
2346
      _fieldName_ = src_iter->first;
 
2347
      if (!fieldMask((int)_fieldName_)) continue; 
 
2348
      DENS_MAT & s(nodalSources[_fieldName_].set_quantity());
 
2349
      myNodalSources[_fieldName_] = DENS_MAT(s.nRows(),s.nCols());
2119
2350
    }
2120
2351
    for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
2121
2352
    {
2162
2393
        }
2163
2394
        
2164
2395
        // assemble
2165
 
        DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
 
2396
        DENS_MAT & s(myNodalSources[_fieldName_].set_quantity());
2166
2397
        _Nmat_ = _fN_.transMat(_fweights_*faceSource);
2167
2398
        for (int i = 0; i < nNodesPerElement_; ++i) 
2168
2399
        {
2169
2400
          int inode = _conn_(i);
2170
2401
          for (int idof = 0; idof < nFieldDOF; ++idof) 
2171
2402
          {
2172
 
            myNodalSources(inode,idof) += _Nmat_(i,idof);
 
2403
            s(inode,idof) += _Nmat_(i,idof);
2173
2404
          }  // end assemble nFieldDOF
2174
2405
        }  // end assemble nNodesPerElement_
2175
2406
      }  // end fset loop
2178
2409
    for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++) {
2179
2410
      _fieldName_ = src_iter->first;
2180
2411
      if (!fieldMask((int)_fieldName_)) continue; 
2181
 
      DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
2182
 
      allsum(communicator_,MPI_IN_PLACE, myNodalSources.ptr(), myNodalSources.size());  
 
2412
      DENS_MAT & s(myNodalSources[_fieldName_].set_quantity());
 
2413
      allsum(communicator_,MPI_IN_PLACE, s.ptr(), s.size());  
 
2414
      DENS_MAT & src(nodalSources[_fieldName_].set_quantity());
 
2415
      src += s;
2183
2416
    }
2184
2417
  }
2185
2418
 
2197
2430
    DENS_MAT elemSource;
2198
2431
    DENS_MAT xCoords(nSD_,nNodesPerElement_);
2199
2432
    DENS_MAT xAtIPs(nSD_,nIPsPerElement_);
 
2433
    FIELDS myNodalSources;
2200
2434
 
2201
 
    if (!(rank_zero(communicator_))) {
2202
 
      // Zero out unmasked nodal sources on all non-main processors.
2203
 
      // This is to avoid counting the previous nodal source values
2204
 
      // multiple times when we aggregate.
2205
 
      for (VOLUME_SOURCE::const_iterator src_iter = sourceFunctions.begin(); 
2206
 
           src_iter != sourceFunctions.end(); src_iter++) {
2207
 
        _fieldName_ = src_iter->first;
2208
 
        int index = (int) _fieldName_;
2209
 
        if (fieldMask(index)) {
2210
 
          DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
2211
 
          myNodalSources.reset(myNodalSources.nRows(), myNodalSources.nCols());
2212
 
        }
2213
 
      }
 
2435
    for (VOLUME_SOURCE::const_iterator src_iter = sourceFunctions.begin(); 
 
2436
         src_iter != sourceFunctions.end(); src_iter++) {
 
2437
      _fieldName_ = src_iter->first;
 
2438
      int index = (int) _fieldName_;
 
2439
      if (! fieldMask(index)) continue;
 
2440
      DENS_MAT & s(nodalSources[_fieldName_].set_quantity());
 
2441
      myNodalSources[_fieldName_] = DENS_MAT(s.nRows(),s.nCols());
2214
2442
    }
2215
 
 
2216
2443
    vector<int> myElems = feMesh_->owned_elts();
2217
 
    for (vector<int>::const_iterator elemsIter = myElems.begin();
 
2444
    for (vector<int>::iterator elemsIter = myElems.begin();
2218
2445
         elemsIter != myElems.end();
2219
2446
         ++elemsIter)
2220
2447
    {
2244
2471
          }
2245
2472
          // assemble
2246
2473
          _Nmat_ = _N_.transMat(_weights_*elemSource);
2247
 
          DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
 
2474
          DENS_MAT & s(myNodalSources[_fieldName_].set_quantity());
2248
2475
          
2249
2476
          for (int i = 0; i < nNodesPerElement_; ++i) {
2250
2477
            int inode = _conn_(i);
2251
2478
            for (int idof = 0; idof < nFieldDOF; ++idof) {
2252
 
              myNodalSources(inode,idof) += _Nmat_(i,idof);
 
2479
              s(inode,idof) += _Nmat_(i,idof);
2253
2480
            }
2254
2481
          }
2255
2482
        }
2260
2487
         src_iter != sourceFunctions.end(); src_iter++) {
2261
2488
      _fieldName_ = src_iter->first;
2262
2489
      int index = (int) _fieldName_;
2263
 
      if (fieldMask(index)) {
2264
 
        DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
2265
 
        allsum(communicator_,MPI_IN_PLACE, myNodalSources.ptr(), myNodalSources.size());  
 
2490
      if (!fieldMask(index)) continue;
 
2491
      DENS_MAT & s(myNodalSources[_fieldName_].set_quantity());
 
2492
      allsum(communicator_,MPI_IN_PLACE, s.ptr(), s.size());  
 
2493
      DENS_MAT & src(nodalSources[_fieldName_].set_quantity());
 
2494
      src += s;
 
2495
    }
 
2496
  }
 
2497
  //-----------------------------------------------------------------
 
2498
  // previously computed nodal sources
 
2499
  //-----------------------------------------------------------------
 
2500
  void FE_Engine::add_sources(const Array<bool> &fieldMask,
 
2501
    const double time,
 
2502
    const FIELDS &sources,
 
2503
    FIELDS &nodalSources) const
 
2504
  {
 
2505
    FIELDS::const_iterator itr;
 
2506
    for (itr=sources.begin(); itr!=sources.end(); itr++) {
 
2507
      _fieldName_ = itr->first;
 
2508
      if (!fieldMask((int)_fieldName_)) continue; 
 
2509
      DENS_MAT & src(nodalSources[_fieldName_].set_quantity());
 
2510
      const DENS_MAT &  s((sources.find(_fieldName_)->second).quantity());
 
2511
      for (int inode = 0; inode < nNodesUnique_; ++inode) {
 
2512
        for (int j = 0; j < src.nCols(); ++j) {
 
2513
          src(inode,j) += s(inode,j);
 
2514
        }
2266
2515
      }
2267
2516
    }
2268
2517
  }
2291
2540
                          // the data member?
2292
2541
 
2293
2542
    // element wise assembly
2294
 
    set< pair <int,int> >::const_iterator iter;
 
2543
    set< pair <int,int> >::iterator iter;
2295
2544
    for (iter = faceSet.begin(); iter != faceSet.end(); iter++) {
2296
2545
      int ielem = iter->first;
2297
2546
      // if this is not our element, do not do calculations
2644
2893
       << feMesh_->num_elements() << " elements";
2645
2894
    print_msg_once(communicator_,ss.str());
2646
2895
#ifdef ATC_VERBOSE
2647
 
    string msg = "element sizes in x:\n";
 
2896
    print_partitions(xmin,xmax,dx);
 
2897
    print_partitions(ymin,ymax,dy);
 
2898
    print_partitions(zmin,zmax,dz);
 
2899
#endif
 
2900
  }
 
2901
  //-----------------------------------------------------------------
 
2902
  void FE_Engine::print_partitions(double xmin, double xmax, Array<double> & dx) const {
 
2903
    stringstream msg;
 
2904
    msg.precision(3);
 
2905
    msg << std::fixed;
 
2906
    msg <<  "\nindex weight fraction location size[A] size[uc]:\n";
2648
2907
    double sum = 0;
2649
2908
    for (int i = 0; i < dx.size(); ++i) { sum += dx(i); }
2650
2909
    double xn= 1.0/sum;
2651
2910
    double xs= xn*(xmax-xmin);
2652
2911
    double xl= xs/(ATC::LammpsInterface::instance()->xlattice());
2653
2912
    double sumn = 0, sums = 0, suml = 0;
 
2913
    double x = xmin;
2654
2914
    for (int i = 0; i < dx.size(); ++i) { 
2655
2915
       double dxn = dx(i)*xn; sumn += dxn;
2656
2916
       double dxs = dx(i)*xs; sums += dxs;
2657
2917
       double dxl = dx(i)*xl; suml += dxl;
2658
 
       msg += "   "+to_string(i+1)
2659
 
             +" "+to_string(dx(i))
2660
 
             +" "+to_string(dxn)
2661
 
             +" "+to_string(dxs)
2662
 
             +" "+to_string(dxl)+"\n";
 
2918
       msg << std::setw(5) << i+1
 
2919
           << std::setw(7) << dx(i)
 
2920
           << std::setw(9) << dxn
 
2921
           << std::setw(9) << x
 
2922
           << std::setw(8) << dxs
 
2923
           << std::setw(9) << dxl << "\n";
 
2924
       x += dxs;
2663
2925
     }
2664
 
     msg += "sum "+to_string(sum)+" "
2665
 
                  +to_string(sumn)+" "
2666
 
                  +to_string(sums)+" "
2667
 
                  +to_string(suml)+"\n";
2668
 
     print_msg_once(communicator_,msg);
2669
 
#endif
 
2926
     msg << "sum  " << setw(7) << sum
 
2927
                    << setw(9) << sumn
 
2928
                    << setw(9) << x
 
2929
                    << setw(8) << sums
 
2930
                    << setw(9) << suml << "\n";
 
2931
     print_msg_once(communicator_,msg.str());
2670
2932
  }
2671
 
 
2672
2933
  //-----------------------------------------------------------------
2673
2934
  // create a uniform rectangular mesh on a rectangular region
2674
2935
  //-----------------------------------------------------------------