From f9d24e3bd05cec52bd0fd2e22d9695e96073153e Mon Sep 17 00:00:00 2001 From: heltai Date: Thu, 25 Apr 2013 17:27:25 +0000 Subject: [PATCH] Merged from trunk to branch_manifold_id git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@29387 0785d39b-7218-0410-832d-ea1e28bc413d --- .../cmake/configure/configure_trilinos.cmake | 26 +- deal.II/include/deal.II/base/logstream.h | 87 +- .../deal.II/matrix_free/fe_evaluation.h | 1400 +++++++++-------- deal.II/source/base/logstream.cc | 34 +- 4 files changed, 780 insertions(+), 767 deletions(-) diff --git a/deal.II/cmake/configure/configure_trilinos.cmake b/deal.II/cmake/configure/configure_trilinos.cmake index 825d304ab3..346d2165fb 100644 --- a/deal.II/cmake/configure/configure_trilinos.cmake +++ b/deal.II/cmake/configure/configure_trilinos.cmake @@ -26,6 +26,13 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var) # SET(${var} TRUE) + # + # Set TRILINOS_DIR to something meaningful if empty + # + IF("${TRILINOS_DIR}" STREQUAL "") + SET(TRILINOS_DIR "") + ENDIF() + # # Check whether all required modules of trilinos are installed: # @@ -34,9 +41,9 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var) ) FOREACH(_module - amesos epetra ifpack aztecoo sacado teuchos + Amesos Epetra Ifpack AztecOO Sacado Teuchos ) - ITEM_MATCHES(_module_found ${_module}$ ${Trilinos_LIBRARIES}) + ITEM_MATCHES(_module_found ${_module} ${Trilinos_PACKAGE_LIST}) IF(_module_found) MESSAGE(STATUS "Found ${_module}") ELSE() @@ -51,12 +58,10 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var) "Missing ${_modules_missing}" ) SET(TRILINOS_ADDITIONAL_ERROR_STRING - "The Trilinos installation found at\n" - " ${TRILINOS_DIR}\n" + "The Trilinos installation (found at \"${TRILINOS_DIR}\")\n" "is missing one or more modules necessary for the deal.II Trilinos interfaces:\n" " ${_modules_missing}\n\n" - "Please re-install Trilinos with the missing Trilinos subpackages - enabled.\n\n" + "Please re-install Trilinos with the missing Trilinos subpackages enabled.\n\n" ) ENDIF() @@ -82,8 +87,7 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var) ) SET(TRILINOS_ADDITIONAL_ERROR_STRING ${TRILINOS_ADDITIONAL_ERROR_STRING} - "The Trilinos installation found at\n" - " ${TRILINOS_DIR}\n" + "The Trilinos installation (found at \"${TRILINOS_DIR}\")\n" "with version ${TRILINOS_VERSION_MAJOR}.${TRILINOS_VERSION_MINOR}.${TRILINOS_VERSION_SUBMINOR} has bugs that make\n" "it incompatible with deal.II. Please use versions before 10.6 or after\n" "10.8.1.\n\n" @@ -103,8 +107,7 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var) ) SET(TRILINOS_ADDITIONAL_ERROR_STRING ${TRILINOS_ADDITIONAL_ERROR_STRING} - "The Trilinos installation found at\n" - " ${TRILINOS_DIR}\n" + "The Trilinos installation (found at \"${TRILINOS_DIR}\")\n" "has to be configured with the same MPI configuration as deal.II, but found:\n" " DEAL_II_WITH_MPI = ${DEAL_II_WITH_MPI}\n" " TRILINOS_WITH_MPI = ${TRILINOS_WITH_MPI}\n" @@ -152,8 +155,7 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var) ) SET(TRILINOS_ADDITIONAL_ERROR_STRING ${TRILINOS_ADDITIONAL_ERROR_STRING} - "The Trilinos installation found at\n" - " ${TRILINOS_DIR}\n" + "The Trilinos installation (found at \"${TRILINOS_DIR}\")\n" "is not compatible with the C++ standard selected for\n" "this compiler. See the deal.II FAQ page for a solution.\n\n" ) diff --git a/deal.II/include/deal.II/base/logstream.h b/deal.II/include/deal.II/base/logstream.h index 5a6166acb6..98617b69a8 100644 --- a/deal.II/include/deal.II/base/logstream.h +++ b/deal.II/include/deal.II/base/logstream.h @@ -531,27 +531,6 @@ private: */ bool at_newline; - /** - * The flags used to modify how - * the stream returned by get_stream() is used. - */ - std::ios_base::fmtflags stream_flags; - - /** - * The minimum width of characters - * used to represent a number. Used by width() to - * change the stream return by get_stream() - */ - std::streamsize stream_width; - - /** - * The maximum width of characters - * used to represent a floating point number. - * Used by precision() to - * change the stream return by get_stream() - */ - std::streamsize stream_precision; - /** * Print head of line. This prints optional time information and the * contents of the prefix stack. @@ -563,43 +542,49 @@ private: * function will return the correct internal ostringstream buffer for * operater<<. */ - std::ostringstream &get_stream() - { - return outstreams.get(); - } + std::ostringstream &get_stream(); /** * We use tbb's thread local storage facility to generate a stringstream * for every thread that sends log messages. */ - Threads::ThreadLocalStorage outstreams; - + Threads::ThreadLocalStorage > outstreams; }; /* ----------------------------- Inline functions and templates ---------------- */ +inline +std::ostringstream & +LogStream::get_stream() +{ + // see if we have already created this stream. if not, do so and + // set the default flags (why we set these flags is lost to + // history, but this is what we need to keep several hundred tests + // from producing different output) + // + // note that in all of this we need not worry about thread-safety + // because we operate on a thread-local object and by definition + // there can only be one access at a time + if (outstreams.get() == 0) + { + outstreams.get().reset (new std::ostringstream); + outstreams.get()->setf(std::ios::showpoint | std::ios::left); + } + + // then return the stream + return *outstreams.get(); +} + + template inline LogStream & LogStream::operator<< (const T &t) { - // save the state of the output stream - std::ostringstream &stream = get_stream(); - - const std::ios::fmtflags old_flags = stream.flags(stream_flags); - const unsigned int old_precision = stream.precision (stream_precision); - const unsigned int old_width = stream.width (stream_width); - - // print to the internal stringstream, using the flags we have just set for - // the stream object: - stream << t; - - // reset output format to whatever it was before - stream.flags (old_flags); - stream.precision(old_precision); - stream.width(old_width); + // print to the internal stringstream + get_stream() << t; return *this; } @@ -610,11 +595,7 @@ inline LogStream & LogStream::operator<< (const double t) { - // save the state of out stream std::ostringstream &stream = get_stream(); - const std::ios::fmtflags old_flags = stream.flags(stream_flags); - const unsigned int old_precision = stream.precision (stream_precision); - const unsigned int old_width = stream.width (stream_width); // we have to make sure that we don't catch NaN's and +-Inf's with the // test, because for these denormals all comparisons are always false. @@ -627,11 +608,6 @@ LogStream::operator<< (const double t) else stream << t*(1.+offset); - // reset output format - stream.flags (old_flags); - stream.precision(old_precision); - stream.width(old_width); - return *this; } @@ -641,11 +617,7 @@ inline LogStream & LogStream::operator<< (const float t) { - // save the state of out stream std::ostringstream &stream = get_stream(); - const std::ios::fmtflags old_flags = stream.flags(stream_flags); - const unsigned int old_precision = stream.precision (stream_precision); - const unsigned int old_width = stream.width (stream_width); // we have to make sure that we don't catch NaN's and +-Inf's with the // test, because for these denormals all comparisons are always false. @@ -658,11 +630,6 @@ LogStream::operator<< (const float t) else stream << t*(1.+offset); - // reset output format - stream.flags (old_flags); - stream.precision(old_precision); - stream.width(old_width); - return *this; } diff --git a/deal.II/include/deal.II/matrix_free/fe_evaluation.h b/deal.II/include/deal.II/matrix_free/fe_evaluation.h index 067488d482..9acb29a336 100644 --- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h +++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h @@ -3789,8 +3789,9 @@ FEEvaluationGeneral namespace internal { // evaluates the given shape data in 1d-3d using the tensor product - // form. does not use the tensor product form and corresponds to a usual - // matrix-matrix product + // form. does not use a particular layout of entries in the matrices + // like the functions below and corresponds to a usual matrix-matrix + // product template inline @@ -3861,6 +3862,724 @@ namespace internal + // This method specializes the general application of tensor-product based + // elements for "symmetric" finite elements, i.e., when the shape functions + // are symmetric about 0.5 and the quadrature points are, too. In that case, + // the 1D shape values read (sorted lexicographically, rows run over 1D + // dofs, columns over quadrature points): + // Q2 --> [ 0.687 0 -0.087 ] + // [ 0.4 1 0.4 ] + // [-0.087 0 0.687 ] + // Q3 --> [ 0.66 0.003 0.002 0.049 ] + // [ 0.521 1.005 -0.01 -0.230 ] + // [-0.230 -0.01 1.005 0.521 ] + // [ 0.049 0.002 0.003 0.66 ] + // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ] + // [ 0.608 1.059 0 0.039 0.176 ] + // [-0.409 -0.113 1 -0.113 -0.409 ] + // [ 0.176 0.039 0 1.059 0.608 ] + // [-0.032 -0.007 0 0.022 0.658 ] + // + // In these matrices, we want to use avoid computations involving zeros and + // ones and in addition use the symmetry in entries to reduce the number of + // read operations. + template + inline + void + apply_tensor_product_values (const VectorizedArray *shape_values, + const VectorizedArray in [], + VectorizedArray out []) + { + AssertIndexRange (direction, dim); + const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, + nn = dof_to_quad ? n_q_points_1d : (fe_degree+1); + const int n_cols = nn / 2; + const int mid = mm / 2; + + const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); + const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); + const int stride = Utilities::fixed_int_power::value; + + for (int i2=0; i2 val0, val1, in0, in1, res0, res1; + if (dof_to_quad == true) + { + val0 = shape_values[col]; + val1 = shape_values[nn-1-col]; + } + else + { + val0 = shape_values[col*n_q_points_1d]; + val1 = shape_values[(col+1)*n_q_points_1d-1]; + } + if (mid > 0) + { + in0 = in[0]; + in1 = in[stride*(mm-1)]; + res0 = val0 * in0; + res1 = val1 * in0; + res0 += val1 * in1; + res1 += val0 * in1; + for (int ind=1; ind(); + if (dof_to_quad == true) + { + if (mm % 2 == 1) + { + val0 = shape_values[mid*n_q_points_1d+col]; + val1 = val0 * in[stride*mid]; + res0 += val1; + res1 += val1; + } + } + else + { + if (mm % 2 == 1 && nn % 2 == 0) + { + val0 = shape_values[col*n_q_points_1d+mid]; + val1 = val0 * in[stride*mid]; + res0 += val1; + res1 += val1; + } + } + if (add == false) + { + out[stride*col] = res0; + out[stride*(nn-1-col)] = res1; + } + else + { + out[stride*col] += res0; + out[stride*(nn-1-col)] += res1; + } + } + if ( dof_to_quad == true && nn%2==1 && mm%2==1 ) + { + if (add==false) + out[stride*n_cols] = in[stride*mid]; + else + out[stride*n_cols] += in[stride*mid]; + } + else if (dof_to_quad == true && nn%2==1) + { + VectorizedArray res0; + VectorizedArray val0 = shape_values[n_cols]; + if (mid > 0) + { + res0 = in[0] + in[stride*(mm-1)]; + res0 *= val0; + for (int ind=1; ind val1 = in[stride*ind] + in[stride*(mm-1-ind)]; + val1 *= val0; + res0 += val1; + } + } + else + res0 = VectorizedArray(); + if (add == false) + out[stride*n_cols] = res0; + else + out[stride*n_cols] += res0; + } + else if (dof_to_quad == false && nn%2 == 1) + { + VectorizedArray res0; + if (mid > 0) + { + VectorizedArray val0 = shape_values[n_cols*n_q_points_1d]; + res0 = in[0] + in[stride*(mm-1)]; + res0 *= val0; + for (int ind=1; ind val1 = in[stride*ind] + in[stride*(mm-1-ind)]; + val1 *= val0; + res0 += val1; + } + if (mm % 2) + res0 += in[stride*mid]; + } + else + res0 = in[0]; + if (add == false) + out[stride*n_cols] = res0; + else + out[stride*n_cols] += res0; + } + + // increment: in regular case, just go to the next point in + // x-direction. If we are at the end of one chunk in x-dir, need to + // jump over to the next layer in z-direction + switch (direction) + { + case 0: + in += mm; + out += nn; + break; + case 1: + case 2: + ++in; + ++out; + break; + default: + Assert (false, ExcNotImplemented()); + } + } + if (direction == 1) + { + in += nn*(mm-1); + out += nn*(nn-1); + } + } + } + + + + // evaluates the given shape data in 1d-3d using the tensor product + // form assuming the symmetries of unit cell shape gradients for + // finite elements in FEEvaluation + + // For the specialized loop used for the gradient computation in + // here, the 1D shape values read (sorted lexicographically, rows + // run over 1D dofs, columns over quadrature points): + // Q2 --> [-2.549 -1 0.549 ] + // [ 3.098 0 -3.098 ] + // [-0.549 1 2.549 ] + // Q3 --> [-4.315 -1.03 0.5 -0.44 ] + // [ 6.07 -1.44 -2.97 2.196 ] + // [-2.196 2.97 1.44 -6.07 ] + // [ 0.44 -0.5 1.03 4.315 ] + // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ] + // [10.111 -2.76 -2.667 2.066 -2.306 ] + // [-5.688 5.773 0 -5.773 5.688 ] + // [ 2.306 -2.066 2.667 2.76 -10.111 ] + // [-0.413 0.353 -0.333 -0.353 0.413 ] + // + // In these matrices, we want to use avoid computations involving + // zeros and ones and in addition use the symmetry in entries to + // reduce the number of read operations. + template + inline + void + apply_tensor_product_gradients (const VectorizedArray *shape_gradients, + const VectorizedArray in [], + VectorizedArray out []) + { + AssertIndexRange (direction, dim); + const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, + nn = dof_to_quad ? n_q_points_1d : (fe_degree+1); + const int n_cols = nn / 2; + const int mid = mm / 2; + + const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); + const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); + const int stride = Utilities::fixed_int_power::value; + + for (int i2=0; i2 val0, val1, in0, in1, res0, res1; + if (dof_to_quad == true) + { + val0 = shape_gradients[col]; + val1 = shape_gradients[nn-1-col]; + } + else + { + val0 = shape_gradients[col*n_q_points_1d]; + val1 = shape_gradients[(nn-col-1)*n_q_points_1d]; + } + if (mid > 0) + { + in0 = in[0]; + in1 = in[stride*(mm-1)]; + res0 = val0 * in0; + res1 = val1 * in0; + res0 -= val1 * in1; + res1 -= val0 * in1; + for (int ind=1; ind(); + if (mm % 2 == 1) + { + if (dof_to_quad == true) + val0 = shape_gradients[mid*n_q_points_1d+col]; + else + val0 = shape_gradients[col*n_q_points_1d+mid]; + val1 = val0 * in[stride*mid]; + res0 += val1; + res1 -= val1; + } + if (add == false) + { + out[stride*col] = res0; + out[stride*(nn-1-col)] = res1; + } + else + { + out[stride*col] += res0; + out[stride*(nn-1-col)] += res1; + } + } + if ( nn%2 == 1 ) + { + VectorizedArray val0, res0; + if (dof_to_quad == true) + val0 = shape_gradients[n_cols]; + else + val0 = shape_gradients[n_cols*n_q_points_1d]; + res0 = in[0] - in[stride*(mm-1)]; + res0 *= val0; + for (int ind=1; ind val1 = in[stride*ind] - in[stride*(mm-1-ind)]; + val1 *= val0; + res0 += val1; + } + if (add == false) + out[stride*n_cols] = res0; + else + out[stride*n_cols] += res0; + } + + // increment: in regular case, just go to the next point in + // x-direction. for y-part in 3D and if we are at the end of one + // chunk in x-dir, need to jump over to the next layer in + // z-direction + switch (direction) + { + case 0: + in += mm; + out += nn; + break; + case 1: + case 2: + ++in; + ++out; + break; + default: + Assert (false, ExcNotImplemented()); + } + } + + if (direction == 1) + { + in += nn * (mm-1); + out += nn * (nn-1); + } + } + } + + + + // evaluates the given shape data in 1d-3d using the tensor product + // form assuming the symmetries of unit cell shape hessians for + // finite elements in FEEvaluation + template + inline + void + apply_tensor_product_hessians (const VectorizedArray *shape_hessians, + const VectorizedArray in [], + VectorizedArray out []) + { + AssertIndexRange (direction, dim); + const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, + nn = dof_to_quad ? n_q_points_1d : (fe_degree+1); + const int n_cols = nn / 2; + const int mid = mm / 2; + + const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); + const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); + const int stride = Utilities::fixed_int_power::value; + + for (int i2=0; i2 val0, val1, in0, in1, res0, res1; + if (dof_to_quad == true) + { + val0 = shape_hessians[col]; + val1 = shape_hessians[nn-1-col]; + } + else + { + val0 = shape_hessians[col*n_q_points_1d]; + val1 = shape_hessians[(col+1)*n_q_points_1d-1]; + } + if (mid > 0) + { + in0 = in[0]; + in1 = in[stride*(mm-1)]; + res0 = val0 * in0; + res1 = val1 * in0; + res0 += val1 * in1; + res1 += val0 * in1; + for (int ind=1; ind(); + if (mm % 2 == 1) + { + if (dof_to_quad == true) + val0 = shape_hessians[mid*n_q_points_1d+col]; + else + val0 = shape_hessians[col*n_q_points_1d+mid]; + val1 = val0 * in[stride*mid]; + res0 += val1; + res1 += val1; + } + if (add == false) + { + out[stride*col] = res0; + out[stride*(nn-1-col)] = res1; + } + else + { + out[stride*col] += res0; + out[stride*(nn-1-col)] += res1; + } + } + if ( nn%2 == 1 ) + { + VectorizedArray val0, res0; + if (dof_to_quad == true) + val0 = shape_hessians[n_cols]; + else + val0 = shape_hessians[n_cols*n_q_points_1d]; + if (mid > 0) + { + res0 = in[0] + in[stride*(mm-1)]; + res0 *= val0; + for (int ind=1; ind val1 = in[stride*ind] + in[stride*(mm-1-ind)]; + val1 *= val0; + res0 += val1; + } + } + else + res0 = VectorizedArray(); + if (mm % 2 == 1) + { + if (dof_to_quad == true) + val0 = shape_hessians[mid*n_q_points_1d+n_cols]; + else + val0 = shape_hessians[n_cols*n_q_points_1d+mid]; + res0 += val0 * in[stride*mid]; + } + if (add == false) + out[stride*n_cols] = res0; + else + out[stride*n_cols] += res0; + } + + // increment: in regular case, just go to the next point in + // x-direction. If we are at the end of one chunk in x-dir, need to + // jump over to the next layer in z-direction + switch (direction) + { + case 0: + in += mm; + out += nn; + break; + case 1: + case 2: + ++in; + ++out; + break; + default: + Assert (false, ExcNotImplemented()); + } + } + if (direction == 1) + { + in += nn*(mm-1); + out += nn*(nn-1); + } + } + } + + + + // evaluates the given shape data in 1d-3d using the tensor product + // form assuming the symmetries of unit cell shape gradients for + // finite elements in FEEvaluationGL + + // This function specializes the application of the tensor product loop for + // Gauss-Lobatto elements which are symmetric about 0.5 just as the general + // class of elements treated by FEEvaluation, have diagonal shape matrices + // for the values and have the following gradient matrices (notice the zeros + // on the diagonal in the interior points, which is due to the construction + // of Legendre polynomials): + // Q2 --> [-3 -1 1 ] + // [ 4 0 -4 ] + // [-1 1 3 ] + // Q3 --> [-6 -1.618 0.618 -1 ] + // [ 8.09 0 -2.236 3.09 ] + // [-3.09 2.236 0 -8.09 ] + // [ 1 -0.618 1.618 6 ] + // Q4 --> [-10 -2.482 0.75 -0.518 1 ] + // [ 13.51 0 -2.673 1.528 -2.82 ] + // [-5.333 3.491 0 -3.491 5.333 ] + // [ 2.82 -1.528 2.673 0 -13.51 ] + // [-1 0.518 -0.75 2.482 10 ] + template + inline + void + apply_tensor_product_gradients_gl (const VectorizedArray *shape_gradients, + const VectorizedArray in [], + VectorizedArray out []) + { + AssertIndexRange (direction, dim); + const int mm = fe_degree+1; + const int nn = fe_degree+1; + const int n_cols = nn / 2; + const int mid = mm / 2; + + const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); + const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); + const int stride = Utilities::fixed_int_power::value; + + for (int i2=0; i2 val0, val1, in0, in1, res0, res1; + if (mid > 0) + { + if (dof_to_quad == true) + { + val0 = shape_gradients[col]; + val1 = shape_gradients[nn-1-col]; + } + else + { + val0 = shape_gradients[col*mm]; + val1 = shape_gradients[(nn-col-1)*mm]; + } + in0 = in[0]; + in1 = in[stride*(mm-1)]; + if (col == 0) + { + if ((mm+dof_to_quad)%2 == 1) + { + res0 = val0 * in0; + res1 = -in0; + res0 += in1; + res1 -= val0 * in1; + } + else + { + res0 = val0 * in0; + res0 -= in1; + res1 = in0; + res1 -= val0 * in1; + } + } + else + { + res0 = val0 * in0; + res1 = val1 * in0; + res0 -= val1 * in1; + res1 -= val0 * in1; + } + for (int ind=1; ind(); + if (mm % 2 == 1) + { + if (dof_to_quad == true) + val0 = shape_gradients[mid*mm+col]; + else + val0 = shape_gradients[col*mm+mid]; + val1 = val0 * in[stride*mid]; + res0 += val1; + res1 -= val1; + } + if (add == false) + { + out[stride*col] = res0; + out[stride*(nn-1-col)] = res1; + } + else + { + out[stride*col] += res0; + out[stride*(nn-1-col)] += res1; + } + } + if ( nn%2 == 1 ) + { + VectorizedArray val0, res0; + if (dof_to_quad == true) + val0 = shape_gradients[n_cols]; + else + val0 = shape_gradients[n_cols*mm]; + if (mid > 0) + { + res0 = in[0] - in[stride*(mm-1)]; + res0 *= val0; + for (int ind=1; ind val1 = in[stride*ind] - in[stride*(mm-1-ind)]; + val1 *= val0; + res0 += val1; + } + } + else + res0 = VectorizedArray(); + if (add == false) + out[stride*n_cols] = res0; + else + out[stride*n_cols] += res0; + } + + // increment: in regular case, just go to the next point in + // x-direction. for y-part in 3D and if we are at the end of one + // chunk in x-dir, need to jump over to the next layer in + // z-direction + switch (direction) + { + case 0: + in += mm; + out += nn; + break; + case 1: + case 2: + ++in; + ++out; + break; + default: + Assert (false, ExcNotImplemented()); + } + } + + if (direction == 1) + { + in += nn * (mm-1); + out += nn * (nn-1); + } + } + } + + + // This performs the evaluation of function values, gradients and Hessians // for tensor-product finite elements. The operation is used for both // FEEvaluationGeneral and FEEvaluation, which provide different functions @@ -4393,8 +5112,6 @@ FEEvaluation -/*----------------- optimized implementation tensor product symmetric case --*/ - template template @@ -4404,196 +5121,9 @@ FEEvaluation ::apply_values (const VectorizedArray in [], VectorizedArray out []) { - AssertIndexRange (direction, dim); - const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, - nn = dof_to_quad ? n_q_points_1d : (fe_degree+1); - const int n_cols = nn / 2; - const int mid = mm / 2; - - const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); - const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); - const int stride = Utilities::fixed_int_power::value; - - // This loop specializes the general application of tensor-product based - // elements for "symmetric" finite elements, i.e., when the shape functions - // are symmetric about 0.5 and the quadrature points are, too. In that case, - // the 1D shape values read (sorted lexicographically, rows run over 1D - // dofs, columns over quadrature points): - // Q2 --> [ 0.687 0 -0.087 ] - // [ 0.4 1 0.4 ] - // [-0.087 0 0.687 ] - // Q3 --> [ 0.66 0.003 0.002 0.049 ] - // [ 0.521 1.005 -0.01 -0.230 ] - // [-0.230 -0.01 1.005 0.521 ] - // [ 0.049 0.002 0.003 0.66 ] - // Q4 --> [ 0.658 0.022 0 -0.007 -0.032 ] - // [ 0.608 1.059 0 0.039 0.176 ] - // [-0.409 -0.113 1 -0.113 -0.409 ] - // [ 0.176 0.039 0 1.059 0.608 ] - // [-0.032 -0.007 0 0.022 0.658 ] - // - // In these matrices, we want to use avoid computations involving zeros and - // ones and in addition use the symmetry in entries to reduce the number of - // read operations. - const VectorizedArray *shape_values = this->data.shape_values.begin(); - for (int i2=0; i2 val0, val1, in0, in1, res0, res1; - if (dof_to_quad == true) - { - val0 = shape_values[col]; - val1 = shape_values[nn-1-col]; - } - else - { - val0 = shape_values[col*n_q_points_1d]; - val1 = shape_values[(col+1)*n_q_points_1d-1]; - } - if (mid > 0) - { - in0 = in[0]; - in1 = in[stride*(mm-1)]; - res0 = val0 * in0; - res1 = val1 * in0; - res0 += val1 * in1; - res1 += val0 * in1; - for (int ind=1; ind(); - if (dof_to_quad == true) - { - if (mm % 2 == 1) - { - val0 = shape_values[mid*n_q_points_1d+col]; - val1 = val0 * in[stride*mid]; - res0 += val1; - res1 += val1; - } - } - else - { - if (mm % 2 == 1 && nn % 2 == 0) - { - val0 = shape_values[col*n_q_points_1d+mid]; - val1 = val0 * in[stride*mid]; - res0 += val1; - res1 += val1; - } - } - if (add == false) - { - out[stride*col] = res0; - out[stride*(nn-1-col)] = res1; - } - else - { - out[stride*col] += res0; - out[stride*(nn-1-col)] += res1; - } - } - if ( dof_to_quad == true && nn%2==1 && mm%2==1 ) - { - if (add==false) - out[stride*n_cols] = in[stride*mid]; - else - out[stride*n_cols] += in[stride*mid]; - } - else if (dof_to_quad == true && nn%2==1) - { - VectorizedArray res0; - VectorizedArray val0 = shape_values[n_cols]; - if (mid > 0) - { - res0 = in[0] + in[stride*(mm-1)]; - res0 *= val0; - for (int ind=1; ind val1 = in[stride*ind] + in[stride*(mm-1-ind)]; - val1 *= val0; - res0 += val1; - } - } - else - res0 = VectorizedArray(); - if (add == false) - out[stride*n_cols] = res0; - else - out[stride*n_cols] += res0; - } - else if (dof_to_quad == false && nn%2 == 1) - { - VectorizedArray res0; - if (mid > 0) - { - VectorizedArray val0 = shape_values[n_cols*n_q_points_1d]; - res0 = in[0] + in[stride*(mm-1)]; - res0 *= val0; - for (int ind=1; ind val1 = in[stride*ind] + in[stride*(mm-1-ind)]; - val1 *= val0; - res0 += val1; - } - if (mm % 2) - res0 += in[stride*mid]; - } - else - res0 = in[0]; - if (add == false) - out[stride*n_cols] = res0; - else - out[stride*n_cols] += res0; - } - - // increment: in regular case, just go to the next point in - // x-direction. If we are at the end of one chunk in x-dir, need to - // jump over to the next layer in z-direction - switch (direction) - { - case 0: - in += mm; - out += nn; - break; - case 1: - case 2: - ++in; - ++out; - break; - default: - Assert (false, ExcNotImplemented()); - } - } - if (direction == 1) - { - in += nn*(mm-1); - out += nn*(nn-1); - } - } + internal::apply_tensor_product_values + (this->data.shape_values.begin(), in, out); } @@ -4607,155 +5137,9 @@ FEEvaluation ::apply_gradients (const VectorizedArray in [], VectorizedArray out []) { - AssertIndexRange (direction, dim); - const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, - nn = dof_to_quad ? n_q_points_1d : (fe_degree+1); - const int n_cols = nn / 2; - const int mid = mm / 2; - - const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); - const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); - const int stride = Utilities::fixed_int_power::value; - - const VectorizedArray *shape_gradients = this->data.shape_gradients.begin(); - for (int i2=0; i2 [-2.549 -1 0.549 ] - // [ 3.098 0 -3.098 ] - // [-0.549 1 2.549 ] - // Q3 --> [-4.315 -1.03 0.5 -0.44 ] - // [ 6.07 -1.44 -2.97 2.196 ] - // [-2.196 2.97 1.44 -6.07 ] - // [ 0.44 -0.5 1.03 4.315 ] - // Q4 --> [-6.316 -1.3 0.333 -0.353 0.413 ] - // [10.111 -2.76 -2.667 2.066 -2.306 ] - // [-5.688 5.773 0 -5.773 5.688 ] - // [ 2.306 -2.066 2.667 2.76 -10.111 ] - // [-0.413 0.353 -0.333 -0.353 0.413 ] - // - // In these matrices, we want to use avoid computations involving - // zeros and ones and in addition use the symmetry in entries to - // reduce the number of read operations. - for (int col=0; col val0, val1, in0, in1, res0, res1; - if (dof_to_quad == true) - { - val0 = shape_gradients[col]; - val1 = shape_gradients[nn-1-col]; - } - else - { - val0 = shape_gradients[col*n_q_points_1d]; - val1 = shape_gradients[(nn-col-1)*n_q_points_1d]; - } - if (mid > 0) - { - in0 = in[0]; - in1 = in[stride*(mm-1)]; - res0 = val0 * in0; - res1 = val1 * in0; - res0 -= val1 * in1; - res1 -= val0 * in1; - for (int ind=1; ind(); - if (mm % 2 == 1) - { - if (dof_to_quad == true) - val0 = shape_gradients[mid*n_q_points_1d+col]; - else - val0 = shape_gradients[col*n_q_points_1d+mid]; - val1 = val0 * in[stride*mid]; - res0 += val1; - res1 -= val1; - } - if (add == false) - { - out[stride*col] = res0; - out[stride*(nn-1-col)] = res1; - } - else - { - out[stride*col] += res0; - out[stride*(nn-1-col)] += res1; - } - } - if ( nn%2 == 1 ) - { - VectorizedArray val0, res0; - if (dof_to_quad == true) - val0 = shape_gradients[n_cols]; - else - val0 = shape_gradients[n_cols*n_q_points_1d]; - res0 = in[0] - in[stride*(mm-1)]; - res0 *= val0; - for (int ind=1; ind val1 = in[stride*ind] - in[stride*(mm-1-ind)]; - val1 *= val0; - res0 += val1; - } - if (add == false) - out[stride*n_cols] = res0; - else - out[stride*n_cols] += res0; - } - - // increment: in regular case, just go to the next point in - // x-direction. for y-part in 3D and if we are at the end of one - // chunk in x-dir, need to jump over to the next layer in - // z-direction - switch (direction) - { - case 0: - in += mm; - out += nn; - break; - case 1: - case 2: - ++in; - ++out; - break; - default: - Assert (false, ExcNotImplemented()); - } - } - - if (direction == 1) - { - in += nn * (mm-1); - out += nn * (nn-1); - } - } + internal::apply_tensor_product_gradients + (this->data.shape_gradients.begin(), in, out); } @@ -4772,146 +5156,9 @@ FEEvaluation ::apply_hessians (const VectorizedArray in [], VectorizedArray out []) { - AssertIndexRange (direction, dim); - const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, - nn = dof_to_quad ? n_q_points_1d : (fe_degree+1); - const int n_cols = nn / 2; - const int mid = mm / 2; - - const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); - const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); - const int stride = Utilities::fixed_int_power::value; - - for (int i2=0; i2 val0, val1, in0, in1, res0, res1; - if (dof_to_quad == true) - { - val0 = this->data.shape_hessians[col]; - val1 = this->data.shape_hessians[nn-1-col]; - } - else - { - val0 = this->data.shape_hessians[col*n_q_points_1d]; - val1 = this->data.shape_hessians[(col+1)*n_q_points_1d-1]; - } - if (mid > 0) - { - in0 = in[0]; - in1 = in[stride*(mm-1)]; - res0 = val0 * in0; - res1 = val1 * in0; - res0 += val1 * in1; - res1 += val0 * in1; - for (int ind=1; inddata.shape_hessians[ind*n_q_points_1d+col]; - val1 = this->data.shape_hessians[ind*n_q_points_1d+nn-1-col]; - } - else - { - val0 = this->data.shape_hessians[col*n_q_points_1d+ind]; - val1 = this->data.shape_hessians[(col+1)*n_q_points_1d-1-ind]; - } - in0 = in[stride*ind]; - in1 = in[stride*(mm-1-ind)]; - res0 += val0 * in0; - res1 += val1 * in0; - res0 += val1 * in1; - res1 += val0 * in1; - } - } - else - res0 = res1 = VectorizedArray(); - if (mm % 2 == 1) - { - if (dof_to_quad == true) - val0 = this->data.shape_hessians[mid*n_q_points_1d+col]; - else - val0 = this->data.shape_hessians[col*n_q_points_1d+mid]; - val1 = val0 * in[stride*mid]; - res0 += val1; - res1 += val1; - } - if (add == false) - { - out[stride*col] = res0; - out[stride*(nn-1-col)] = res1; - } - else - { - out[stride*col] += res0; - out[stride*(nn-1-col)] += res1; - } - } - if ( nn%2 == 1 ) - { - VectorizedArray val0, res0; - if (dof_to_quad == true) - val0 = this->data.shape_hessians[n_cols]; - else - val0 = this->data.shape_hessians[n_cols*n_q_points_1d]; - if (mid > 0) - { - res0 = in[0] + in[stride*(mm-1)]; - res0 *= val0; - for (int ind=1; inddata.shape_hessians[ind*n_q_points_1d+n_cols]; - else - val0 = this->data.shape_hessians[n_cols*n_q_points_1d+ind]; - VectorizedArray val1 = in[stride*ind] + in[stride*(mm-1-ind)]; - val1 *= val0; - res0 += val1; - } - } - else - res0 = VectorizedArray(); - if (mm % 2 == 1) - { - if (dof_to_quad == true) - val0 = this->data.shape_hessians[mid*n_q_points_1d+n_cols]; - else - val0 = this->data.shape_hessians[n_cols*n_q_points_1d+mid]; - res0 += val0 * in[stride*mid]; - } - if (add == false) - out[stride*n_cols] = res0; - else - out[stride*n_cols] += res0; - } - - // increment: in regular case, just go to the next point in - // x-direction. If we are at the end of one chunk in x-dir, need to - // jump over to the next layer in z-direction - switch (direction) - { - case 0: - in += mm; - out += nn; - break; - case 1: - case 2: - ++in; - ++out; - break; - default: - Assert (false, ExcNotImplemented()); - } - } - if (direction == 1) - { - in += nn*(mm-1); - out += nn*(nn-1); - } - } + internal::apply_tensor_product_hessians + (this->data.shape_hessians.begin(), in, out); } @@ -5152,193 +5399,14 @@ FEEvaluationGL ::apply_gradients (const VectorizedArray in [], VectorizedArray out []) { - AssertIndexRange (direction, dim); - const int mm = fe_degree+1; - const int nn = fe_degree+1; - const int n_cols = nn / 2; - const int mid = mm / 2; - - const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); - const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); - const int stride = Utilities::fixed_int_power::value; - - // This loop specializes the application of the tensor product loop for - // Gauss-Lobatto elements which are symmetric about 0.5 just as the general - // class of elements treated by FEEvaluation, have diagonal shape matrices - // for the values and have the following gradient matrices (notice the zeros - // on the diagonal in the interior points, which is due to the construction - // of Legendre polynomials): - // Q2 --> [-3 -1 1 ] - // [ 4 0 -4 ] - // [-1 1 3 ] - // Q3 --> [-6 -1.618 0.618 -1 ] - // [ 8.09 0 -2.236 3.09 ] - // [-3.09 2.236 0 -8.09 ] - // [ 1 -0.618 1.618 6 ] - // Q4 --> [-10 -2.482 0.75 -0.518 1 ] - // [ 13.51 0 -2.673 1.528 -2.82 ] - // [-5.333 3.491 0 -3.491 5.333 ] - // [ 2.82 -1.528 2.673 0 -13.51 ] - // [-1 0.518 -0.75 2.482 10 ] - for (int i2=0; i2 val0, val1, in0, in1, res0, res1; - if (mid > 0) - { - if (dof_to_quad == true) - { - val0 = this->data.shape_gradients[col]; - val1 = this->data.shape_gradients[nn-1-col]; - } - else - { - val0 = this->data.shape_gradients[col*mm]; - val1 = this->data.shape_gradients[(nn-col-1)*mm]; - } - in0 = in[0]; - in1 = in[stride*(mm-1)]; - if (col == 0) - { - if ((mm+dof_to_quad)%2 == 1) - { - res0 = val0 * in0; - res1 = -in0; - res0 += in1; - res1 -= val0 * in1; - } - else - { - res0 = val0 * in0; - res0 -= in1; - res1 = in0; - res1 -= val0 * in1; - } - } - else - { - res0 = val0 * in0; - res1 = val1 * in0; - res0 -= val1 * in1; - res1 -= val0 * in1; - } - for (int ind=1; inddata.shape_gradients[ind*mm+col]; - val1 = this->data.shape_gradients[ind*mm+nn-1-col]; - } - else - { - val0 = this->data.shape_gradients[col*mm+ind]; - val1 = this->data.shape_gradients[(nn-col-1)*mm+ind]; - } - - // at inner points, the gradient is zero for ind==col - in0 = in[stride*ind]; - in1 = in[stride*(mm-1-ind)]; - if (ind == col) - { - res1 += val1 * in0; - res0 -= val1 * in1; - } - else - { - res0 += val0 * in0; - res1 += val1 * in0; - res0 -= val1 * in1; - res1 -= val0 * in1; - } - } - } - else - res0 = res1 = VectorizedArray(); - if (mm % 2 == 1) - { - if (dof_to_quad == true) - val0 = this->data.shape_gradients[mid*mm+col]; - else - val0 = this->data.shape_gradients[col*mm+mid]; - val1 = val0 * in[stride*mid]; - res0 += val1; - res1 -= val1; - } - if (add == false) - { - out[stride*col] = res0; - out[stride*(nn-1-col)] = res1; - } - else - { - out[stride*col] += res0; - out[stride*(nn-1-col)] += res1; - } - } - if ( nn%2 == 1 ) - { - VectorizedArray val0, res0; - if (dof_to_quad == true) - val0 = this->data.shape_gradients[n_cols]; - else - val0 = this->data.shape_gradients[n_cols*mm]; - if (mid > 0) - { - res0 = in[0] - in[stride*(mm-1)]; - res0 *= val0; - for (int ind=1; inddata.shape_gradients[ind*mm+n_cols]; - else - val0 = this->data.shape_gradients[n_cols*mm+ind]; - VectorizedArray val1 = in[stride*ind] - in[stride*(mm-1-ind)]; - val1 *= val0; - res0 += val1; - } - } - else - res0 = VectorizedArray(); - if (add == false) - out[stride*n_cols] = res0; - else - out[stride*n_cols] += res0; - } - - // increment: in regular case, just go to the next point in - // x-direction. for y-part in 3D and if we are at the end of one - // chunk in x-dir, need to jump over to the next layer in - // z-direction - switch (direction) - { - case 0: - in += mm; - out += nn; - break; - case 1: - case 2: - ++in; - ++out; - break; - default: - Assert (false, ExcNotImplemented()); - } - } + internal::apply_tensor_product_gradients_gl + (this->data.shape_gradients.begin(), in, out); +} - if (direction == 1) - { - in += nn * (mm-1); - out += nn * (nn-1); - } - } #endif // ifndef DOXYGEN -} - DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/base/logstream.cc b/deal.II/source/base/logstream.cc index 7594dc2750..588f1af187 100644 --- a/deal.II/source/base/logstream.cc +++ b/deal.II/source/base/logstream.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009, 2010, 2011, 2012 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009, 2010, 2011, 2012, 2013 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -57,10 +57,7 @@ LogStream::LogStream() float_threshold(0.), offset(0), old_cerr(0), - at_newline(true), - stream_flags(std::ios::showpoint | std::ios::left), - stream_width(std::cout.width()), - stream_precision(std::cout.precision()) + at_newline(true) { get_prefixes().push("DEAL:"); @@ -134,10 +131,6 @@ LogStream::operator<< (std::ostream& (*p) (std::ostream &)) { std::ostringstream &stream = get_stream(); - // save the state of out stream - std::ios::fmtflags old_flags = stream.flags(stream_flags); - unsigned int old_precision = stream.precision (stream_precision); - unsigned int old_width = stream.width (stream_width); // Print to the internal stringstream: stream << p; @@ -170,11 +163,6 @@ LogStream::operator<< (std::ostream& (*p) (std::ostream &)) stream.str(""); } - // reset output format - stream.flags (old_flags); - stream.precision(old_precision); - stream.width(old_width); - return *this; } @@ -260,33 +248,21 @@ void LogStream::pop () std::ios::fmtflags LogStream::flags(const std::ios::fmtflags f) { - Threads::Mutex::ScopedLock lock(log_lock); - - std::ios::fmtflags tmp = stream_flags; - stream_flags = f; - return tmp; + return get_stream().flags (f); } std::streamsize LogStream::precision (const std::streamsize prec) { - Threads::Mutex::ScopedLock lock(log_lock); - - std::streamsize tmp = stream_precision; - stream_precision = prec; - return tmp; + return get_stream().precision (prec); } std::streamsize LogStream::width (const std::streamsize wide) { - Threads::Mutex::ScopedLock lock(log_lock); - - std::streamsize tmp = stream_width; - stream_width = wide; - return tmp; + return get_stream().width (wide); } -- 2.39.5