Re: [deal.II] ILU preconditioner and mesh-scale dependency?

[Jaekwang Kim]

Thanks for sharing insight! 


On Thursday, March 14, 2019 at 12:18:30 PM UTC-5, Wolfgang Bangerth wrote:
>
> On 3/13/19 8:18 PM, Jaekwang Kim wrote: 
> > 
> > 
> > I wonder now ... 
> > 
> > 1. why smaller scale mesh only results such error, 
>
> You have two terms in your matrix, one from the advection term and one 
> from the reaction term. If you scale the mesh by a factor of c, the 
> height of the shape functions used in the reaction term does not change, 
> but the gradient of shape functions used in the advection term is scaled 
> by 1/c. As a consequence, the balance of the two terms that together 
> make up the matrix changes and you apparently get into a situation where 
> the resulting matrix is ill-conditioned or indeed not invertible. 
>
>
> > 2. If possible, is there any choice of better preconditionner for my 
> > system other than SparseILU to go around this problem? 
>
> To make sure that the system you have assembled is correct, I usually 
> first check with the SparseDirectUMFPACK solver. Once I know that the 
> system is correct, I think about iterative solvers and preconditioners. 
> In your case, the system you have is advection dominated, for which it 
> is notoriously difficult to construct good preconditioners. 
>
> Best 
>   WB 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] ILU preconditioner and mesh-scale dependency?

[Jaekwang Kim]

Hi 

I had a trouble in solving scalar transport equation using DG elements. 

For example, I was solving for \lambda in the following equation where *u *is 
also a numerical solution

[image: Screen Shot 2019-03-13 at 9.08.17 PM.png]


The thing is I do not get any problem when my mesh file is scale of O(1) 
and I have tested my code works fine using the Method of Manufactured 
solution. 

Yet, when I tried to solve the governing equation at small scale 
mesh...(I.e. all (x,y) point in the original mesh-file was reduced by 
(0.01x,0.01y) 

I receive error like this 

An error occurred in line <129> of file 

 
in function

void dealii::SparseILU::initialize(const 
SparseMatrix &, const dealii::SparseILU::AdditionalData &) 
[number = double, somenumber = double]

The violated condition was: 

luval[ia[k]] != 0

Additional information: 

While computing the ILU decomposition, the algorithm found a zero pivot 
on the diagonal of row 34. This must stop the ILU algorithm because it 
means that the matrix for which you try to compute a decomposition is 
singular.


The problem occurs when the red line of the code is executed...

 template 

void StokesProblem::solve_transport ()

{

std::cout << "   -solve transport begins"<< std::endl;

SolverControl   solver_control (std::pow(10,6), 
system_rhs.block(2).l2_norm() * pow(10,-4));



unsigned int restart = 500;

SolverGMRES< Vector >::AdditionalData 
gmres_additional_data(restart+2);

SolverGMRES< Vector > solver(solver_control, 
gmres_additional_data);

  

//make preconditioner

SparseILU::AdditionalData additional_data(0,500); // (0 , 
additional diagonal terms)

std::cout << "B" << std::endl;

SparseILU preconditioner;



*preconditioner.initialize (system_matrix.block(2,2), 
additional_data);*



std::cout << "A" << std::endl;

solver.solve (system_matrix.block(2,2), solution.block(2), 
system_rhs.block(2), preconditioner);

// constraints.distribute (solution);

}




I wonder now ...

1. why smaller scale mesh only results such error,

2. If possible, is there any choice of better preconditionner for my system 
other than SparseILU to go around this problem? 


Thanks,

Regards, 

Jaekwang Kim 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Install error when compile with Trilinos

[JAEKWANG KIM]

Thanks
I complied it successfully

On Sun, Feb 24, 2019 at 10:55 PM Bruno Turcksin 
wrote:

> Jaekwang,
>
> This is the problem:
>
> Le dim. 24 févr. 2019 à 22:10, Jaekwang Kim  a
> écrit :
>
> > -- Include
> /Users/jaekwangkim/Program/dealii-9.0.1/cmake/configure/configure_2_trilinos.cmake
> >
> > CMake Warning at /usr/local/lib/cmake/Trilinos/TrilinosConfig.cmake:156
> (MESSAGE):
> >
> >   Component "Pike" NOT found.
> >
> > Call Stack (most recent call first):
> >
> >   cmake/macros/macro_find_package.cmake:27 (_FIND_PACKAGE)
> >
> >   cmake/modules/FindTRILINOS.cmake:41 (FIND_PACKAGE)
> >
> >   cmake/macros/macro_find_package.cmake:27 (_FIND_PACKAGE)
> >
> >   cmake/configure/configure_2_trilinos.cmake:22 (FIND_PACKAGE)
> >
> >
>  
> /Users/jaekwangkim/Program/dealii/CMakeFiles/CMakeTmp/evaluate_expression.tmp:1
> (FEATURE_TRILINOS_FIND_EXTERNAL)
> >
> >   cmake/macros/macro_evaluate_expression.cmake:30 (INCLUDE)
> >
> >   cmake/macros/macro_configure_feature.cmake:237 (EVALUATE_EXPRESSION)
> >
> >   cmake/configure/configure_2_trilinos.cmake:219 (CONFIGURE_FEATURE)
> >
> >   cmake/macros/macro_verbose_include.cmake:19 (INCLUDE)
> >
> >   CMakeLists.txt:124 (VERBOSE_INCLUDE)
> >
> >
> >
> > CMake Error at /usr/local/lib/cmake/Trilinos/TrilinosConfig.cmake:219
> (include):
> >
> >   include could not find load file:
> >
> >
> > /usr/local/lib/cmake/Trilinos/../Pike/PikeConfig.cmake
>
> Basically what happens is that Trilinos says that Pike is installed
> but when deal.II is trying to find the configuration file for Pike, it
> can't find it. These kind of errors used to happen for the old
> versions of Trilinos but I haven't seen that in a while. Try to
> install a newer version of Trilinos and if you don't need Pike, don't
> enable it.
>
> Best,
>
> Bruno
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/r0pKGTxIjrI/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Install error when compile with Trilinos

[Jaekwang Kim]

I resummarized screen messages that seems to be a problem 



...~~

*-- Performing Test DEAL_II_HAVE_AVX*

*-- Performing Test DEAL_II_HAVE_AVX - Failed*

*-- Performing Test DEAL_II_HAVE_AVX512*

*-- Performing Test DEAL_II_HAVE_AVX512 - Failed*

 

*DEAL_II_COMPILER_DEFAULTS_TO_CXX11_OR_NEWER - Failed*

*DEAL_II_HAVE_COMPLEX_OPERATOR_OVERLOADS - Failed*

 

 

*-- Performing Test gfortran_LIBRARY*

*-- Performing Test gfortran_LIBRARY - Failed*

*-- Performing Test gcc_ext.10.5_LIBRARY*

*-- Performing Test gcc_ext.10.5_LIBRARY - Failed*

*-- Performing Test gcc_LIBRARY*

*-- Performing Test gcc_LIBRARY - Failed*

*-- Performing Test quadmath_LIBRARY*

*-- Performing Test quadmath_LIBRARY - Failed*

*-- Include /Users/jaekwangkim/Program/dealii-*

*-- Performing Test DEAL_II_HAVE_MT_POSIX_BARRIERS*

*-- Performing Test DEAL_II_HAVE_MT_POSIX_BARRIERS - Failed*

-- Performing Test DEAL_II_HAVE_FLAG_Wno_parentheses

-- Performing Test DEAL_II_HAVE_FLAG_Wno_parentheses - Success

-- DEAL_II_WITH_THREADS successfully set up with bundled packages.

-- 

-- Performing Test DEAL_II_HAS_AUTO_PTR

-- Performing Test DEAL_II_HAS_AUTO_PTR - Failed

 

CMake Error at /usr/local/lib/cmake/Trilinos/TrilinosConfig.cmake:219 
(include):

  include could not find load file:

 

   /usr/local/lib/cmake/Trilinos/../Pike/PikeConfig.cmake

 

-- Found EPETRA_CONFIG_H

-- Found SACADO_CONFIG_H

-- Found SACADO_CMATH_HPP

-- Performing Test TRILINOS_SUPPORTS_CPP11

-- Performing Test TRILINOS_SUPPORTS_CPP11 - Success

-- Performing Test TRILINOS_HAS_C99_TR1_WORKAROUND

-- Performing Test TRILINOS_HAS_C99_TR1_WORKAROUND - Success

-- Found SACADO_TRAD_HPP

-- Performing Test TRILINOS_CXX_SUPPORTS_SACADO_COMPLEX_RAD

-- Performing Test TRILINOS_CXX_SUPPORTS_SACADO_COMPLEX_RAD - Failed

-- Found TRILINOS_LIBRARY_trilinoscouplings

-- Found TRILINOS_LIBRARY_msqutil

-- Found TRILINOS_LIBRARY_mesquite

-- Found TRILINOS_LIBRARY_ctrilinos

-- Found TRILINOS_LIBRARY_rol

-- Found TRILINOS_LIBRARY_moochothyra

-- Found TRILINOS_LIBRARY_moocho

-- Found TRILINOS_LIBRARY_rythmos

-- Found TRILINOS_LIBRARY_muelu-adapters

-- Found TRILINOS_LIBRARY_muelu-interface

-- Found TRILINOS_LIBRARY_muelu

-- Found TRILINOS_LIBRARY_moertel

-- Found TRILINOS_LIBRARY_phalanx

-- Found TRILINOS_LIBRARY_intrepid

-- Found TRILINOS_LIBRARY_ifpack2-adapters

-- Found TRILINOS_LIBRARY_ifpack2

-- Found TRILINOS_LIBRARY_anasazitpetra

-- Found TRILINOS_LIBRARY_ModeLaplace

-- Found TRILINOS_LIBRARY_anasaziepetra

-- Found TRILINOS_LIBRARY_anasazi

-- Found TRILINOS_LIBRARY_komplex

-- Found TRILINOS_LIBRARY_suplib_cpp

-- Found TRILINOS_LIBRARY_suplib_c

-- Found TRILINOS_LIBRARY_suplib

-- Found TRILINOS_LIBRARY_supes

-- Found TRILINOS_LIBRARY_aprepro_lib

-- Found TRILINOS_LIBRARY_chaco

-- Found TRILINOS_LIBRARY_io_info_lib

-- Found TRILINOS_LIBRARY_Ionit

-- Found TRILINOS_LIBRARY_Iotr

-- Found TRILINOS_LIBRARY_Iohb

-- Found TRILINOS_LIBRARY_Iogn

-- Found TRILINOS_LIBRARY_Iovs

-- Found TRILINOS_LIBRARY_Iopg

-- Found TRILINOS_LIBRARY_Ioss

-- Found TRILINOS_LIBRARY_amesos2

-- Found TRILINOS_LIBRARY_shylu

-- Found TRILINOS_LIBRARY_belostpetra

-- Found TRILINOS_LIBRARY_belosepetra

-- Found TRILINOS_LIBRARY_belos

-- Found TRILINOS_LIBRARY_ml

-- Found TRILINOS_LIBRARY_ifpack

-- Found TRILINOS_LIBRARY_zoltan2

-- Found TRILINOS_LIBRARY_pamgen_extras

-- Found TRILINOS_LIBRARY_pamgen

-- Found TRILINOS_LIBRARY_amesos

-- Found TRILINOS_LIBRARY_galeri-xpetra

-- Found TRILINOS_LIBRARY_galeri-epetra

-- Found TRILINOS_LIBRARY_aztecoo

-- Found TRILINOS_LIBRARY_dpliris

-- Found TRILINOS_LIBRARY_isorropia

-- Found TRILINOS_LIBRARY_optipack

-- Found TRILINOS_LIBRARY_xpetra-sup

-- Found TRILINOS_LIBRARY_xpetra

-- Found TRILINOS_LIBRARY_thyratpetra

-- Found TRILINOS_LIBRARY_thyraepetraext

-- Found TRILINOS_LIBRARY_thyraepetra

-- Found TRILINOS_LIBRARY_thyracore

-- Found TRILINOS_LIBRARY_domi

-- Found TRILINOS_LIBRARY_epetraext

-- Found TRILINOS_LIBRARY_trilinosss

-- Found TRILINOS_LIBRARY_tpetraext

-- Found TRILINOS_LIBRARY_tpetrainout

-- Found TRILINOS_LIBRARY_tpetra

-- Found TRILINOS_LIBRARY_kokkostsqr

-- Found TRILINOS_LIBRARY_tpetraclassiclinalg

-- Found TRILINOS_LIBRARY_tpetraclassicnodeapi

-- Found TRILINOS_LIBRARY_tpetraclassic

-- Found TRILINOS_LIBRARY_triutils

-- Found TRILINOS_LIBRARY_globipack

-- Found TRILINOS_LIBRARY_shards

-- Found TRILINOS_LIBRARY_zoltan

-- Found TRILINOS_LIBRARY_simpi

-- Found TRILINOS_LIBRARY_epetra

-- Found TRILINOS_LIBRARY_minitensor

-- Found TRILINOS_LIBRARY_sacado

-- Found TRILINOS_LIBRARY_rtop

-- Found TRILINOS_LIBRARY_kokkoskernels

-- Found TRILINOS_LIBRARY_teuchoskokkoscomm

-- Found TRILINOS_LIBRARY_teuchoskokkoscompat

-- Found TRILINOS_LIBRARY_teuchosremainder

-- Found TRILINOS_LIBRARY_teuchosnumerics

-- Found TRILINOS_LIBRARY_teuchoscomm

-- Found TRILINOS_LIBRARY_teuchosparameterlist

-- Found TRILINOS_LIBRARY_teuchoscore

-- Found TRILINOS_LIBRARY_kokkosalgorithms

[deal.II] Install error when compile with Trilinos

[Jaekwang Kim]

Hi, all. 

I have used deal.ii for year, and recently I need to compiled deal.ii with 
Trilinos 

I installed Trilinos on my local enabling only serial (non MPI) 

However, after the installation, when I tried to reinstall with deal.ii 
enabling Trilinos packages, 

I got error so that I cannot 'make install' 

My configuration is as follow and I also attach my 'CMakeError.log' file. 

I would appreciate any idea to resolve present issue

Thanks!

Jaekwang 

###

#

#  deal.II configuration:

#CMAKE_BUILD_TYPE:   DebugRelease

#BUILD_SHARED_LIBS:  ON

#CMAKE_INSTALL_PREFIX:   /Users/jaekwangkim/Program/dealii

#CMAKE_SOURCE_DIR:   /Users/jaekwangkim/Program/dealii-9.0.1

#(version 9.0.1)

#CMAKE_BINARY_DIR:   /Users/jaekwangkim/Program/dealii

#CMAKE_CXX_COMPILER: AppleClang 10.0.0.10001044 on platform 
Darwin x86_64

#
/Library/Developer/CommandLineTools/usr/bin/c++

#

#  Configured Features (DEAL_II_ALLOW_BUNDLED = ON, 
DEAL_II_ALLOW_AUTODETECTION = ON):

#  ( DEAL_II_WITH_64BIT_INDICES = OFF )

#  ( DEAL_II_WITH_ADOLC = OFF )

#  ( DEAL_II_WITH_ARPACK = OFF )

#  ( DEAL_II_WITH_ASSIMP = OFF )

#DEAL_II_WITH_BOOST set up with external dependencies

#  ( DEAL_II_WITH_CUDA = OFF )

#DEAL_II_WITH_CXX14 = ON

#DEAL_II_WITH_CXX17 = ON

#  ( DEAL_II_WITH_GMSH = OFF )

#  ( DEAL_II_WITH_GSL = OFF )

#  ( DEAL_II_WITH_HDF5 = OFF )

#DEAL_II_WITH_LAPACK set up with external dependencies

#  ( DEAL_II_WITH_METIS = OFF )

#  ( DEAL_II_WITH_MPI = OFF )

#DEAL_II_WITH_MUPARSER set up with bundled packages

#  ( DEAL_II_WITH_NANOFLANN = OFF )

#  ( DEAL_II_WITH_NETCDF = OFF )

#  ( DEAL_II_WITH_OPENCASCADE = OFF )

#  ( DEAL_II_WITH_P4EST = OFF )

#  ( DEAL_II_WITH_PETSC = OFF )

#  ( DEAL_II_WITH_SCALAPACK = OFF )

#  ( DEAL_II_WITH_SLEPC = OFF )

#  ( DEAL_II_WITH_SUNDIALS = OFF )

#DEAL_II_WITH_THREADS set up with bundled packages

#DEAL_II_WITH_TRILINOS set up with external dependencies

#DEAL_II_WITH_UMFPACK set up with external dependencies

#DEAL_II_WITH_ZLIB set up with external dependencies

#

#  Component configuration:

#  ( DEAL_II_COMPONENT_DOCUMENTATION = OFF )

#DEAL_II_COMPONENT_EXAMPLES

#  ( DEAL_II_COMPONENT_PACKAGE = OFF )

#  ( DEAL_II_COMPONENT_PYTHON_BINDINGS = OFF )

#

#  Detailed information (compiler flags, feature configuration) can be 
found in detailed.log

#

#  Run  $ make info  to print a help message with a list of top level 
targets

#

###

*-- Configuring incomplete, errors occurred!*




-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Performing C++ SOURCE FILE Test DEAL_II_HAVE_FLAG_Wl__as_needed failed with the 
following output:
Change Dir: /Users/jaekwangkim/Program/dealii/CMakeFiles/CMakeTmp

Run Build Command:"/usr/bin/make" "cmTC_76208/fast"
/Library/Developer/CommandLineTools/usr/bin/make -f 
CMakeFiles/cmTC_76208.dir/build.make CMakeFiles/cmTC_76208.dir/build
Building CXX object CMakeFiles/cmTC_76208.dir/src.cxx.o
/Library/Developer/CommandLineTools/usr/bin/c++
-DDEAL_II_HAVE_FLAG_Wl__as_needed -isysroot 
/Library/Developer/CommandLineTools/SDKs/MacOSX10.14.sdk   -o 
CMakeFiles/cmTC_76208.dir/src.cxx.o -c 
/Users/jaekwangkim/Program/dealii/CMakeFiles/CMakeTmp/src.cxx
Linking CXX executable cmTC_76208
/usr/local/bin/cmake -E cmake_link_script CMakeFiles/cmTC_76208.dir/link.txt 
--verbose=1
/Library/Developer/CommandLineTools/usr/bin/c++   
-DDEAL_II_HAVE_FLAG_Wl__as_needed -isysroot 
/Library/Developer/CommandLineTools/SDKs/MacOSX10.14.sdk 
-Wl,-search_paths_first -Wl,-headerpad_max_install_names   
CMakeFiles/cmTC_76208.dir/src.cxx.o  -o cmTC_76208 -Wl,--as-needed 
ld: unknown option: --as-needed
clang: error: linker command failed with exit code 1 (use -v to see invocation)
make[1]: *** [cmTC_76208] Error 1
make: *** [cmTC_76208/fast] Error 2

Source file was:
int main() { return 0; }
Performing C++ SOURCE FILE Test DEAL_II_HAVE_FLAG_Wno_placement_new failed with 
the following output:
Change Dir: /Users/jaekwangkim/Program/dealii/CMakeFiles/CMakeTmp

Run Build Command:"/usr/bin/make" "cmTC_6a621/fast"
/Library/Developer/CommandLineTools/usr/bin/make -f 
CMakeFiles/cmTC_6a621.dir/build.make CMakeFiles/cmTC_6a621.dir/build
Building CXX object CMakeFiles/cmTC_6a621.dir/src.cxx.o
/Library/Developer/CommandLineTools/usr/bin/c++
-DDEAL_II_HAVE_FLAG_Wno_placement_new -isysroot 

[deal.II] Re: Discrete Galekin with periodic boundary condition

[Jaekwang Kim]

No, you were right, 
I couldn't get the correct results without 
'Triangulation::add_periodicity()'. 

The current my local environment does not configure deal.ii with MPI and 
Trilinos as the example step-45, which shows how this function is used. 
But I guess I can still use the "  'Triangulation::add_periodicity() " 
 function 

For example, 

 template 

void AdvectionProblem::create_grid ()

{

double pi = 3.14159265359;

const Point<2> end (4*pi,1);

const Point<2> start (0.0,0.0);



std::vector< unsigned int > repetitions(2);

repetitions[0]=5;

repetitions[1]=1;

 

GridGenerator::subdivided_hyper_rectangle 
(triangulation,repetitions, start, end, false);



for (typename Triangulation::active_cell_iterator

 cell=triangulation.begin_active();

 cell!=triangulation.end(); ++cell)

{

cell->set_refine_flag(RefinementCase::cut_x);



for (unsigned int f=0; f::faces_per_cell; ++f)

{

// You need this line to ensure that this doesn't affect 
anything other than bondary faces!

if (cell->face(f)->at_boundary() == true)

{

cell->face(f)->set_boundary_id(0);



// Boundary faces

static const double tol = 1e-12;

if ( std::abs(cell->face(f)->center()[0]-0.0)< tol )

{

cell->face(f)->set_boundary_id(5);

std::cout << " read boundary 5 " << std::endl;

}

if ( std::abs(cell->face(f)->center()[0]-4*pi)< tol )

{

cell->face(f)->set_boundary_id(6);

std::cout << " read boundary 6 " << std::endl;

}

if ( std::abs(cell->face(f)->center()[1]-0.0)< tol )

{

cell->face(f)->set_boundary_id(2);

//std::cout << " read boundary 2 " << std::endl;

}

if ( std::abs(cell->face(f)->center()[1]-1.0)< tol )

{

cell->face(f)->set_boundary_id(4);

//std::cout << " read boundary 4 " << std::endl;

}

}

}

}



//original argument

std::vector::cell_iterator>> pairs;



 

GridTools::collect_periodic_faces(dof_handler,

  /*b_id1*/ 6,

  /*b_id2*/ 5,

  /*direction*/ 0,

  pairs);

triangulation.add_periodicity(pairs);



triangulation.refine_global (3);



}





However, I got error when I use add_periodicity function which says that 


no viable conversion from


'vector, false>::CellAccessor>>>' to 'const


vector>>>'


triangulation.add_periodicity(pairs);


^





candidate constructor not viable: no known conversion from


'std::vector::cell_iterator>


>' (aka 'vector, false> >


> >') to 'const


std::__1::vector > >,


std::__1::allocator > > > > &' for 1st argument


vector(const vector& __x);


^


Do you have any idea how to fix the problem ? 


Thanks!


Best, 

Jaekwang 



On Friday, November 9, 2018 at 7:46:02 PM UTC-6, Daniel Arndt wrote:
>
> Jaekwang,
>
> It seems that I don't need to add_periodicity to the triangulation unless 
>> it is parallel distributed triangulation? 
>>
> I would be heavily surprised if you get correct results without calling 
> Triangulation::add_periodicity() if you are using discontinuous elements 
> and a serial Triangulation.
> Otherwise, you are not providing any kind of information to these objects 
> that you want to treat specific boundaries as periodic ones. 
>  
>
>> By the way, when I 'add_periodicity', I have compiler error on the red 
>>> line of the code. 
>>>
>>> Can you give me a clue to resolve this error? 
>>>
>>>   //original argument
>>>
>>> std::vector>>
>>> PeriodicFacePair::cell_iterator>> 
>>> pairs;
>>>
>>> 
>>>
>>> GridTools::collect_periodic_faces(dof_handler,
>>>
>>>   /*b_id1*/ 6,
>>>
>>>   /*b_id2*/ 5,
>>>
>>>   /*direction*/ 0,
>>>
>>>   pairs);
>>>
>>> 
>>>
>>> triangulation.add_periodicity(pairs);
>>>
>> You need to call GridTools::collect_periodic_faces with a parameter of 
> type 
>
> std::vector Triangulation::cell_iterator>> pairs;
>
> if you want to use the resulting object for 
> Triangulation::add_periodicity().
>  
> Best,
> Daniel
>

-- 

[deal.II] Re: Discrete Galekin with periodic boundary condition

[Jaekwang Kim]

You were right. I couldn't get the correct results without 
 Triangulation::add_periodicity(). 

So, it seems that I should compile deal.ii with MPI,P4EST,Trilinos to use 
that function(Triangulation::add_periodicity()) since my triangulation 
should be parallel distributed...

Right? 

Best, 
Jaekwang 
 

On Friday, November 9, 2018 at 7:46:02 PM UTC-6, Daniel Arndt wrote:
>
> Jaekwang,
>
> It seems that I don't need to add_periodicity to the triangulation unless 
>> it is parallel distributed triangulation? 
>>
> I would be heavily surprised if you get correct results without calling 
> Triangulation::add_periodicity() if you are using discontinuous elements 
> and a serial Triangulation.
> Otherwise, you are not providing any kind of information to these objects 
> that you want to treat specific boundaries as periodic ones. 
>  
>
>> By the way, when I 'add_periodicity', I have compiler error on the red 
>>> line of the code. 
>>>
>>> Can you give me a clue to resolve this error? 
>>>
>>>   //original argument
>>>
>>> std::vector>>
>>> PeriodicFacePair::cell_iterator>> 
>>> pairs;
>>>
>>> 
>>>
>>> GridTools::collect_periodic_faces(dof_handler,
>>>
>>>   /*b_id1*/ 6,
>>>
>>>   /*b_id2*/ 5,
>>>
>>>   /*direction*/ 0,
>>>
>>>   pairs);
>>>
>>> 
>>>
>>> triangulation.add_periodicity(pairs);
>>>
>> You need to call GridTools::collect_periodic_faces with a parameter of 
> type 
>
> std::vector Triangulation::cell_iterator>> pairs;
>
> if you want to use the resulting object for 
> Triangulation::add_periodicity().
>  
> Best,
> Daniel
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Discrete Galekin with periodic boundary condition

[Jaekwang Kim]

It seems that I don't need to add_periodicity to the triangulation unless 
it is parallel distributed triangulation? 


On Friday, November 9, 2018 at 9:45:53 AM UTC-6, Jaekwang Kim wrote:
>
> Thanks! Now it becomes clear to me. 
>
> By the way, when I 'add_periodicity', I have compiler error on the red 
> line of the code. 
>
> Can you give me a clue to resolve this error? 
>
>   //original argument
>
> std::vector
> PeriodicFacePair::cell_iterator>> pairs;
>
> 
>
> GridTools::collect_periodic_faces(dof_handler,
>
>   /*b_id1*/ 6,
>
>   /*b_id2*/ 5,
>
>   /*direction*/ 0,
>
>   pairs);
>
> 
>
> triangulation.add_periodicity(pairs);
>
> Error messages says that 
>
>
> *no viable conversion from*
>
> *  'vector *dealii::internal::DoFHandler::Iterators
> *  2>, false>::CellAccessor**>>>' to '**const*
>
>   *vector**>>>'*
>
> triangulation.add_periodicity(pairs);
>
> *  ^*
>
>
> also says that 
>
>
>   candidate constructor not viable: no known conversion from
>
>   'std::vector DoFHandler<2>::cell_iterator>
>
>   >' (aka 
> 'vector, 
> false> >
>
>   > >') to 'const
>
>   
> std::__1::vector
>   2> > >,
>
>   
> std::__1::allocator
>   2> > > > > &' for 1st argument
>
> vector(const vector& __x);
>
> *^*
>
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Discrete Galekin with periodic boundary condition

[Jaekwang Kim]

Thanks! Now it becomes clear to me. 

By the way, when I 'add_periodicity', I have compiler error on the red line 
of the code. 

Can you give me a clue to resolve this error? 

  //original argument

std::vector::cell_iterator>> pairs;



GridTools::collect_periodic_faces(dof_handler,

  /*b_id1*/ 6,

  /*b_id2*/ 5,

  /*direction*/ 0,

  pairs);



triangulation.add_periodicity(pairs);

Error messages says that 


*no viable conversion from*

*  'vector, false>::CellAccessor**>>>' to '**const*

  *vector**>>>'*

triangulation.add_periodicity(pairs);

*  ^*


also says that 


  candidate constructor not viable: no known conversion from

  'std::vector::cell_iterator>

  >' (aka 
'vector, 
false> >

  > >') to 'const

  
std::__1::vector > >,

  
std::__1::allocator > > > > &' for 1st argument

vector(const vector& __x);

*^*



-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Discrete Galekin with periodic boundary condition

[Jaekwang Kim]

Thanks, for the reply..

Hi am using following function to construct periodicity when I set up dof 



DoFTools::make_hanging_node_constraints (dof_handler, constraints);

DoFTools::make_periodicity_constraints(dof_handler,6,5 //Boundary ID 6 is 
outlet and 5 is inlet 

   ,0,

constraints

);

 


Question 1. 

I assume that this does same thing as Triangulation.add_periodicity () + 
Gridtools::collect_periodic_faces(). 

Isn't it? 


Question 2. 

What I don't know now is when I call 'integrate_boundary_term' function in 
the mesh worker 

how can I assess neighbor cell information?

For example, when I assemble for left end cell - the influx should be 
described by right-end cell at previous time solution solution at previous 
time.

To summary, I don't know how to call neighbor (described by periodic) cell 
solution within the 'integrate_boundary term' function



Thanks !! 

Jaekwang 


   template 

void AdvectionProblem::integrate_boundary_term (DoFInfo ,

 CellInfo )

{ 


FullMatrix _matrix = dinfo.matrix(0).matrix;

Vector _vector = dinfo.vector(0).block(0);



const FEValuesBase _v = info.fe_values();

const std::vector  = fe_v.get_JxW_values ();

const std::vector >  = 
fe_v.get_all_normal_vectors ();



const AdvectionField advection_field;

std::vector > advection_values 
(fe_v.n_quadrature_points);

advection_field.value_list 
(fe_v.get_quadrature_points(),advection_values);












On Thursday, November 8, 2018 at 2:16:11 AM UTC-6, Daniel Arndt wrote:
>
> Jaekwang,
>
> In general, using periodic boundary conditions with MeshWorker should be 
> pretty easy.
> All you would have to do is calling Triangulation.add_periodicity() 
> 
> with a parameter constructed by GridTools::collect_periodic_faces() 
> 
> .
> Doing so will let MeshWorker treat all periodic faces as inner faces.
>
> The 1D case is less well tested though. If you encounter additional 
> difficulties feel free to ask.
>
> Best,
> Daniel
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Discrete Galekin with periodic boundary condition

[Jaekwang Kim]

Thanks for the reply! 
I will take a look

Jaekwang 

On Wednesday, November 7, 2018 at 10:05:10 PM UTC-6, Praveen C wrote:
>
> You have to do some extra work to use periodic bc with MeshWorker. I have 
> a Euler DG code here which can do that
>
> https://github.com/cpraveen/dflo/tree/master/src_mpi
>
> The periodic bc is setup in these files (not my code, authors are 
> mentioned in the files)
>
> DealiiExtensions.cc
> DealiiExtensions.h 
> <https://github.com/cpraveen/dflo/blob/master/src_mpi/DealiiExtensions.h>
>
> Then see how these functions are used in claw.cc
>
> This code is somewhat old now and I have not run it for quite some time. 
> So it may or may not compile/run with latest deal.II, you have to try it 
> out.
>
> I have a simpler DG code for linear advection without MeshWorker here
>
> https://bitbucket.org/cpraveen/deal_ii/src/master/dg/1d_scalar_legendre/
>
> Best
> praveen
>
> On 08-Nov-2018, at 8:31 AM, Jaekwang Kim > 
> wrote:
>
>
>
> I am trying to solve advection equation using DG method with Meshworker. 
>
> Yet, I have a problem in implementing periodic boundary condition. 
>
> For example, I consider a 1D advection equation with periodic boundary 
> condition
>
>
>
> with source term f=-sin(x), I assume that I should get cos(x) as a steady 
> state solution, and I have confirmed this set up for the toy problem using 
> Finite Difference.  
>
>
>
>
> I derived my weak form as follow...
>
> 
>
>
>
> I am assembling equation (11) every time step,
>
>
> Mass matrix for n+1 step solution 
>
>
>  local_matrix(i,j) += fe_v.shape_value(i,point) *
>
>   fe_v.shape_value(j,point) *
>
>   JxW[point];
>
>
>
>
>
> and right hand side is 
>
> (cell integration) 
>
> local_vector(i) += (fe_v.shape_value(i,point)
>
> *(sol_phi[point] + dt * source_values[point])
>
> +
>
> dt*(advection_values[point]
>
> *fe_v.shape_grad (i,point))
>
> *sol_phi[point]
>
> ) *JxW[point];
>
>
>
> (Face integration) 
>
>
>  for (unsigned int i=0; i
> {
>
> if(advection_dot_normal >0) //outflux
>
> {local_vector(i) -= dt * advection_dot_normal *
>
> 
> (fe_v.shape_value(i,point)-fe_v_neighbor.shape_value(i,point)) *
>
> my_sol_phi[point] *
>
> JxW[point];
>
> }
>
> else //influx
>
> {
>
> local_vector(i) -= dt * advection_dot_normal *
>
>
> (fe_v.shape_value(i,point)-fe_v_neighbor.shape_value(i,point)) *
>
>neighbor_sol_phi[point] *
>
>JxW[point];
>
> }
>
> }
>
>
> (Boundary Integration)  <-think this is may be the problem 
>
>
>
> for (unsigned int i=0; i
> {   //sign changed here
>
> if(advection_dot_normal >0) //outflux
>
> { local_vector(i)-= dt*advection_dot_normal *
>
> fe_v.shape_value(i,point) *
>
> sol_phi[point] *
>
> JxW[point];
>
> }
>
> else
>
> {
>
> //influx this is problematic...it should be 
> connected to last end cell (but..how?)
>
>   local_vector(i)-= dt*advection_dot_normal *
>
> fe_v.shape_value(i,point) *
>
> *sol_phi[point] * *it should be 
> connected to last end cell, but the boundary integrate function cannot 
> call information from reference cell...  (but..how?)
>
> JxW[point];
>
> }
>
> }
>
>
> As I am solving this x in [0,4pi] I should see two period of cosine wave 
> after enough time step, but it is not..
>
>
>
> 
>
>
> Is there anyone who have the idea to use Meshworker for 
>
>
> *1D advection equation with periodic boundary condition?*
>
>
>
> Thanks in advance 
>
>
> Jaekwang 
>
>
>
> -- 
> T

[deal.II] Stokes flow with Periodic boundary condition

[Jaekwang Kim]

Hi, all. 

I'd like ask a question implementing periodic boundary condition. 

Modifying step-22 tutorial code, I want to solve a stokes flow at Couette 
configuration

I plan to solve it at (x,y) in [0,5] x [0,1], rectangular domain 

with boundary condition 

\vec{u} (x,1)= (1,0) (moving dragging plate)
\vec{u} (x,0)= (0,0) (stationary bottom plate) 

and periodic boundary condition at both inlet and outlet 

\vec{u}(0,y)=\vec{u}(5,y)



I think in deal.ii, this can be implemented as 

{

  constraints.clear ();


  FEValuesExtractors::Vector velocities(0);

  DoFTools::make_hanging_node_constraints (dof_handler,

   constraints);

  

  // PERIODIC boundary condition



  std::vector::cell_iterator> > make_pair;



  GridTools::collect_periodic_faces(dof_handler,1,3

  /*direction*/,0

  ,make_pair);



  VectorTools::interpolate_boundary_values  (dof_handler,

4,

BoundaryValues(),

constraints,


fe.component_mask(velocities));



  VectorTools::interpolate_boundary_values (dof_handler,

2,

ZeroFunction(dim+1),

constraints,


fe.component_mask(velocities));



}



However, as I run code, I get error message that 

An error occurred in line <3751> of file 
 in 
function

void dealii::GridTools::collect_periodic_faces(const MeshType &, const 
types::boundary_id, const types::boundary_id, const int, 
std::vector > &, const 
Tensor<1, MeshType::space_dimension> &, const FullMatrix &) 
[MeshType = dealii::DoFHandler<2, 2>]

The violated condition was: 

pairs1.size() == pairs2.size()

Additional information: 

Unmatched faces on periodic boundaries


Stacktrace:

---

when it tried to compile  following part 

* GridTools::collect_periodic_faces(dof_handler,1,3*

*  /*direction*/,0*

*  ,make_pair);*



I would appreciate any comments or thoughts for possible reason of the 
problem.

I attach minimalized problem that has same problem. 



Thanks



-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include  //For the periodic boundary condition

#include 
#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 
#include 

#include 

#include 

#include 
#include 
#include 

namespace Step22
{
  using namespace dealii;


  template 
  struct InnerPreconditioner;

  template <>
  struct InnerPreconditioner<2>
  {
typedef SparseDirectUMFPACK type;
  };

  template <>
  struct InnerPreconditioner<3>
  {
typedef SparseILU type;
  };



  template 
  class StokesProblem
  {
  public:
StokesProblem (const unsigned int degree);
void run ();

  private:
void define_mesh ();
void setup_dofs ();
void assemble_system ();
void solve ();
void output_results ();

const unsigned int   degree;

Triangulation   triangulation;
FESystemfe;
DoFHandler  dof_handler;

ConstraintMatrix constraints;

BlockSparsityPattern  sparsity_pattern;
BlockSparseMatrix system_matrix;

BlockVector solution;
BlockVector system_rhs;

std_cxx11::shared_ptr::type> A_preconditioner;
  };



  template 
  class BoundaryValues : public Function
  {
  public:
BoundaryValues () : Function(dim+1) {}

virtual double value (const Point   ,
  const unsigned int  component = 0) const;

virtual void vector_value (const Point ,
   Vector   ) const;
  };


  template 
  double
  BoundaryValues::value (const Point  ,
  const unsigned int component) const
  {
Assert (component < this->n_components,
ExcIndexRange (component, 0, this->n_components));

if (component == 0)
return 1;
return 0;
  }


  template 
  void
  BoundaryValues::vector_value (const Point ,
 

Re: [deal.II] interpolating boundary values from a previously found FE solution

[Jaekwang Kim]

Dr. Bangerth

Is there any example tutorial related to the first approach? 


*"You can build a ConstraintMatrix object that sets boundary DoFs of V to 
the corresponding values of U. This should be easy to do." *

On Friday, March 30, 2018 at 8:15:23 AM UTC-5, Wolfgang Bangerth wrote:
>
> On 03/29/2018 11:18 PM, Mustafa Aggul wrote: 
> > 
> > suppose I assemble a system and solve it. Say my solution is U. In the 
> next 
> > step, I want to solve a different problem for V. My desire is to make V 
> > exactly same as U on the boundaries. 
> > 
>
> You can build a ConstraintMatrix object that sets boundary DoFs of V to 
> the 
> corresponding values of U. This should be easy to do. 
>
> Alternatively, you could solve for a function D so that V=U+D. D should 
> then 
> simply have zero boundary values. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Adding rows on system matrix connecting to dof

[Jaekwang Kim]

Hi, all 

I have a question on using deal.ii for solving unknown boundary condition 
problem. 
I have thought the most minimized related to my question. 

where I consider a navies slip boundary condition, and I have tried to 
solve this with finite difference.
For more demonstration please see following screenshot. 




Finite Difference Formulation might be...





where the last raw of matrix is added to show boundary condition at lower 
plate.

 

and I think problem is generally well posed so that I can solve this with 







Now, I am trying to solve same problem using finite element method...





The description of same boundary condition in weak form will look like 

 

 has on the boundary



for (unsigned int face=0; faceface(face)->boundary_id() == 2)

{  for (unsigned int q=0; q

[deal.II] Re: Use of 'VectorTools::interpolate_boundary_values' using numerical solution.

[Jaekwang Kim]

Thanks for giving me an idea for this!
I will take a look approaches you just mentioned. 

Jaekwang

On Wednesday, April 11, 2018 at 12:24:26 AM UTC-5, Jie Cheng wrote:
>
> Hi Jaekwang
>
> One way I can think of is to construct your ConstraintMatrix for Part B 
> manually. Iterate over your triangulation, find all the dofs on the target 
> boundaries, then use ConstraintMatrix::add_line and 
> ConstraintMatrix::set_inhomogeinity to set the BC.
>
> Another method might be putting your calc_tangential_velocity function 
> inside New_VelBndValues. I assume this function needs to know the solution, 
> triangulation etc. Then you would have to create some references or dealii 
> SmartPointers in this class as well.
>
> I hope this helps.
>
> Jie
>
> On Monday, April 9, 2018 at 6:32:28 PM UTC-4, Jaekwang Kim wrote:
>>
>> Hi, all 
>>
>>
>> I am always thank you for all guys here in this community. 
>>
>>
>> I have a question on the use of  VectorTools::interpolate_boundary_values 
>> functions
>>
>>
>> Typically we use this function in when we set-up system as in the form 
>>
>>
>>
>> VectorTools::interpolate_boundary_values (dof_handler,
>>
>> 1,
>>
>>*Original**_VelBndValues(),*
>>
>>constraints,
>>
>>fe.component_mask(velocities));
>>
>>
>>
>> The red lined functions are usually declared out side of main class. 
>>
>>
>> One problem I am interested in now is to (Part A) solve stokes flow with 
>> no slip boundary condition on every boundary one time (as step-22), -
>>
>>
>> and then to (Part B) solve same system again with slightly different 
>> boundary condition comes from the straight problem
>>
>>
>> If I finish Part A, i have numerical solution of (u,p) linked with dot 
>> handler. 
>>
>>
>> In (Part B), I post process solutions part A, e.g. I calculated stress 
>> value and I calculated tangential velocity of a boundary using Navier slip 
>> condition, 
>>
>>
>> So, I need to use 
>>
>>
>>
>> VectorTools::interpolate_boundary_values (dof_handler,
>>
>> 1,
>>
>>*New_VelBndValues(),*
>>
>>constraints,
>>
>>fe.component_mask(velocities));
>>
>>
>> such that 
>>
>> 
>>
>>  template
>>
>>  double
>>
>>  New_VelBndValues::value (constPoint  ,
>>
>>  constunsignedintcomponent) const
>>
>>  {
>>
>>  Assert (component < this->n_components,
>>
>>  ExcIndexRange (component, 0, this->n_components));
>>
>>  
>>
>>  if(component == 0)
>>
>>  *return calc_tangential_velocity (); (Error Here!!)*
>>
>>  else
>>
>>  return0;
>>
>>  }
>>
>>
>>
>>
>> where '*calc_tangential_velocity ();' *is declared in main class, since 
>> I should have information of dot_handler and solution vectors.
>>
>>
>>   template
>>
>>   classStokesProblem
>>
>>   {
>>
>>   public:
>>
>> StokesProblem (constunsignedintdegree);
>>
>> voidrun ();
>>
>>   private:
>>
>> 
>>
>>   ……
>>
>>  
>>
>> *double calc_tangential_velocity ();*
>>
>>  
>>
>> BlockVector solution;
>>
>>  
>>
>> };
>>
>> 
>>
>>   
>>
>> Obviously, this method does not work, I mean I cannot compile the code 
>> and it shows  'use of undeclared identifier 'calc_tangential_velocity' on 
>> the line  *(Error Here!!)*
>>
>> I would appreciate any one who comment it. It seems that complier has to 
>> know what is that before. 
>>
>>
>> Will this be possible in a certain way? How can I  go around this problem 
>> ? 
>>
>>
>> Thanks, 
>>
>>
>> Jaekwang 
>>
>>
>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Use of 'VectorTools::interpolate_boundary_values' using numerical solution.

[Jaekwang Kim]

Hi, all 


I am always thank you for all guys here in this community. 


I have a question on the use of  VectorTools::interpolate_boundary_values 
functions


Typically we use this function in when we set-up system as in the form 



VectorTools::interpolate_boundary_values (dof_handler,

1,

   *Original**_VelBndValues(),*

   constraints,

   fe.component_mask(velocities));



The red lined functions are usually declared out side of main class. 


One problem I am interested in now is to (Part A) solve stokes flow with no 
slip boundary condition on every boundary one time (as step-22), -


and then to (Part B) solve same system again with slightly different 
boundary condition comes from the straight problem


If I finish Part A, i have numerical solution of (u,p) linked with dot 
handler. 


In (Part B), I post process solutions part A, e.g. I calculated stress 
value and I calculated tangential velocity of a boundary using Navier slip 
condition, 


So, I need to use 



VectorTools::interpolate_boundary_values (dof_handler,

1,

   *New_VelBndValues(),*

   constraints,

   fe.component_mask(velocities));


such that 



 template

 double

 New_VelBndValues::value (constPoint  ,

 constunsignedintcomponent) const

 {

 Assert (component < this->n_components,

 ExcIndexRange (component, 0, this->n_components));

 

 if(component == 0)

 *return calc_tangential_velocity (); (Error Here!!)*

 else

 return0;

 }




where '*calc_tangential_velocity ();' *is declared in main class, since I 
should have information of dot_handler and solution vectors.


  template

  classStokesProblem

  {

  public:

StokesProblem (constunsignedintdegree);

voidrun ();

  private:



  ……

 

*double calc_tangential_velocity ();*

 

BlockVector solution;

 

};



  

Obviously, this method does not work, I mean I cannot compile the code and 
it shows  'use of undeclared identifier 'calc_tangential_velocity' on the 
line  *(Error Here!!)*

I would appreciate any one who comment it. It seems that complier has to 
know what is that before. 


Will this be possible in a certain way? How can I  go around this problem ? 


Thanks, 


Jaekwang 



-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: use of "meshworker" for FESystem

[JAEKWANG KIM]

Thanks for the attetion.
I overcame this problem, though I don’t recognize what was solution now.
It turns out to be possible…

Thank you for your time

Jaekwang



On Sun, Mar 25, 2018 at 3:13 AM, Daniel Arndt <
daniel.ar...@iwr.uni-heidelberg.de> wrote:

> Jaekwang,
>
> do you have a minimal compilign example that shows the problem?
> Did you call dof_handler.distribute_dofs(fe)?
>
> Best,
> Daniel
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/
> forum/dealii?hl=en
> ---
> You received this message because you are subscribed to a topic in the
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/
> topic/dealii/zDT4ItovbJ4/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] use of "meshworker" for FESystem

[Jaekwang Kim]

Hi, all.

I am curious about whether i can use  MeshWorker: class for FESystem.

I am trying to solve system, 

that I need to solve 'stokes flow' and then use velocity solution to solve 
other scalar transport equation using DG method 

that looks like 


I solve for \lambda, and u is known from stokes solver. h(x,y) and f(x,y) 
is also known.  


So I declared my FESystem as, 


  template 

StokesProblem::StokesProblem (const unsigned int degree)

:

degree (degree),

triangulation (Triangulation::allow_anisotropic_smoothing),

fe (FE_Q(degree+1), dim,

FE_Q(degree), 1,

FE_DGQ(degree+1), 1),

dof_handler (triangulation)

{}


Afte solving stokes system just as in step-22, 

I need mesh-worker after that 

however when I try to reinitialize info_box, I receive errors like 

 template 

void StokesProblem::setup_mesh_worker ()

{

std::cout << "Setting up mesh worker ...\n";



MeshWorker::IntegrationInfoBox info_box;

const unsigned int n_gauss_points = dof_handler.get_fe().degree+1; // 
I haven't thought about this.



//Since quadratures for cells, boundary, and interior faces can be 
selected independently,

//We have to hand over this three times.

info_box.initialize_gauss_quadrature(n_gauss_points,

 n_gauss_points,

 n_gauss_points);



info_box.initialize_update_flags();

UpdateFlags update_flags = update_quadrature_points |

update_values|

update_gradients;

info_box.add_update_flags(update_flags, true, true, true, true);



const MappingQ mapping (degree);




info_box.initialize (fe , mapping); //Note that  ' FESystem fe'

}


I am receiving error saying that


*An error occurred in line <1198> of file 
 in 
function*

*const FiniteElement ::DoFHandler<2, 2>::get_fe() 
const [dim = 2, spacedim = 2]*

*The violated condition was: *

*selected_fe!=0*

*The name and call sequence of the exception was:*

*ExcMessage("You are trying to access the DoFHandler's FiniteElement 
object before it has been initialized.")*

*Additional Information: *

*You are trying to access the DoFHandler's FiniteElement object before it 
has been initialized.*

I would appreciate very much anyone who helps me...

Thanks, 

Jaekwang 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Relation between Solution Error Behavior and Polynomial Approximation Degree

[Jaekwang Kim]

No, I typically, calculate values itself and use MATLAB to draw plots. 

Thanks,

Jaekwang 

On Monday, November 27, 2017 at 3:24:46 PM UTC-6, seven wrote:
>
> Hello Jaekwang,
>
> I am trying to generate some log-log plots, and wondering if you used the 
> functions in deal.ii to generate the figure. If not, what did you use?
>
> Thanks,
> Jiaqi
>
> 在 2016年9月29日星期四 UTC-4上午11:41:48，Jaekwang Kim写道：
>>
>>
>> <https://lh3.googleusercontent.com/-dkfwJvgPEmc/V-00Y-c4tzI/A5s/pb0Z-IBRe-An_S39dSDzxWg1k5GJlEbKQCLcB/s1600/step7plotting.jpg>
>>
>> Hi all, I have question on error behavior of FEM. 
>>
>> I thought that the order of error is O(h^p) where h is a mesh-size and p 
>> is polynomial degree we use in approximation. 
>>
>> So, I thought that if I plot an error with number of mesh in log-log 
>> scale, than the graph will show -p slope. 
>> However, I the error behaves little bit different from my expectation.
>>
>> For example, I use a step7 tutorial program (which solves Helmholtz 
>> decomposition and compares the FEM solution with exact solution.) 
>>
>> The error curve showed more steep slope whenever I increase polynomial 
>> degree approximation however, the slope is not (-p). 
>> I reached slope (-3) when I used fifth-degree polynomial approximation... 
>>   
>> You can check this behavior in attached picture. 
>>
>> Until now, I have considered, 
>>
>> 1. Mapping(From reference cell to real cell) degree (which is originally 
>> set to 1 but I used higher mapping) 
>> 2. Instead of Qgauss quadrature, I am using QgaussLobatto Quadrature for 
>> any integration over cells. 
>> 3. Shape function , again I tried to use QgaussLobatto node point for 
>> this) 
>>
>> is there any suggestion that I need to fix more? 
>> or my first prediction that the slope will show '-p' or error will just 
>> behave O(h^p) was wrong?
>>
>> I am always thank you for all guys!
>>
>> Jaekwang Kim  
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Distribute constraints on PETSc wrapper

[Jaekwang Kim]

Hi, all 

I was practicing using Petsc library through deal.ii wrappers, 

While, step 17 gives me clear explanation, I am having trouble in 
distributing hanging node constraints on Petsc system matrix, 

I was modifying step-9 hyperbolic equation solvers with PETSc matrix, but 
when I compile my code  I run into error message like 

*[0]PETSC ERROR: - Error Message 
*

*[0]PETSC ERROR: Floating point exception!*

*[0]PETSC ERROR: Inserting nan at matrix entry (0,0)!*

*[0]PETSC ERROR: 
*

*[0]PETSC ERROR: Petsc Release Version 3.3.0, Patch 6, Mon Feb 11 12:26:34 
CST 2013 *

*[0]PETSC ERROR: See docs/changes/index.html for recent updates.*

*[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.*

*[0]PETSC ERROR: See docs/index.html for manual pages.*

*[0]PETSC ERROR: 
*

*[0]PETSC ERROR: ./advection on a dbg named golubh2 by jk12 Thu Sep 21 
13:05:46 2017*

*[0]PETSC ERROR: Libraries linked from 
/usr/local/apps/petsc/3.3-p6/intel/mvapich2-1.6-intel/dbg/lib*

*[0]PETSC ERROR: Configure run at Wed Mar 13 14:51:28 2013*

*[0]PETSC ERROR: Configure options 
--prefix=/usr/local/apps/petsc/3.3-p6/intel/mvapich2-1.6-intel/dbg 
--with-cc=mpicc --with-fc=mpif90 --with-cxx=mpicxx 
--with-shared-libraries=0 --with-mpi=1 --known-mpi-shared=1 
--with-blas-lapack-lib="[/usr/local/intel-11.1/mkl/lib/em64t/libmkl_intel_lp64.a,/usr/local/intel-11.1/mkl/lib/em64t/libmkl_sequential.a,/usr/local/intel-11.1/mkl/lib/em64t/libmkl_core.a]"*

*[0]PETSC ERROR: 
*

*[0]PETSC ERROR: MatSetValues() line 1013 in 
/usr/local/apps/petsc/build/petsc-3.3-p6/src/mat/interface/matrix.c*

*ERROR: Uncaught exception in MPI_InitFinalize on proc 0. Skipping 
MPI_Finalize() to avoid a deadlock.*

I have checked where I got stuck in my code and found that they are 
generated when I distribute constraints to (yellow highlighted lines). 

I attached how I made constraints matrix at setup system. 

I would appreciate if any one can suggests where I got stuck... 

Thanks

Jaekwang 





 template 

void AdvectionProblem::assemble_system_petsc ()

{

pcout << "assemble system A" << std::endl;



const MappingQ mapping (degree);

QGauss   quadrature_formula(3);

QGauss face_quadrature_formula(3);



FEValues fe_values (mapping, fe, quadrature_formula,

 update_values   | update_gradients |

 update_quadrature_points | 
update_JxW_values);



FEFaceValues fe_face_values (mapping, fe, 
face_quadrature_formula,

  update_values| 
update_normal_vectors | update_gradients |

  update_quadrature_points  | 
update_JxW_values);



const AdvectionField advection_field; //Defined by Tensor 
Function

const RightHandSide  right_hand_side; //Defined by normal 
Vector Function

const BoundaryValues boundary_values; //Defined by normal 
Vector Function



const unsigned int   dofs_per_cell = fe.dofs_per_cell;

const unsigned int   n_q_points= quadrature_formula.size();

const unsigned int   n_face_q_points = 
fe_face_values.get_quadrature().size();



FullMatrix   cell_matrix (dofs_per_cell, dofs_per_cell);

Vector   cell_rhs (dofs_per_cell);



std::vector local_dof_indices 
(dofs_per_cell);



std::vector rhs_values (n_q_points);

std::vector > advection_directions (n_q_points);

std::vector face_boundary_values (n_face_q_points);

std::vector > face_advection_directions 
(n_face_q_points);



double tau ;



pcout << "assemble system B" << std::endl;



typename DoFHandler::active_cell_iterator

cell = dof_handler.begin_active(),

endc = dof_handler.end();

for (; cell!=endc; ++cell)

if (cell->subdomain_id() == this_mpi_process)

{



pcout << "assemble system C" << std::endl;



cell_matrix = 0;

cell_rhs = 0;



tau =  0.5 * cell->diameter () * cell->diameter ();



advection_field.value_list 
(fe_values.get_quadrature_points(),

advection_directions);

right_hand_side.value_list 
(fe_values.get_quadrature_points(),

rhs_values);





// pcout << "assemble system D " << 

[deal.II] Step-9, iterator over the cell

[Jaekwang Kim]

Hi, all. 
I have a questions on an algorithm in step-9 that does some iteration over 
the cell with its neighbors. 

I tried to read explanation carefully, but some part is unclear yet. 

In step-9, one class estimates a gradient of solution in each cell by 
calculating estimated value from its neighbor cells. 

Special care might be needed when our domain consist of different level of 
refinement in order to include all neighbor cells, and I think this parts 
from step-9 does this... but I am not understanding this.. 

I made some example, in the simple mesh above, 


*case 1: *
when we iterate over on cell-4 in my pic, it has four faces and each face 
has its neighbor, (2,3,5,6) 

for cell, 2 and 3, it has same refinement level and I assume that 2 and 3 
become automatically active. 

and I wonder what happens when it comes to 5 and 6. where the neighbor cell 
has one coarse level? 

The explanation in step 9 tutorial tells us that.. 

*If we are not in 1d, we collect all neighbor children `behind' the 
subfaces of the current face*
for (unsigned int subface_no=0; subface_non_children(); ++subface_no)
active_neighbors.push_back (
std::get<0>(*cell)->neighbor_child_on_subface(face_no,subface_no));
}
}


and I think I don't understand the blue line properlyand especially the 
word subfaces means.

*case 2: *

when we iterative over on cell-5 in my pic, I am assuming that the cell 2,4 
and 7 are automatically activated neighbor cell... 

and we don't need to anything special.. 

I would appreciate if someone can help me. 

Thanks!!


Jaekwang 



 

template 

void

GradientEstimation::estimate_cell (const 
SynchronousIterators > ,

   EstimateScratchData 
  _data,

   const EstimateCopyData &)

{

Tensor<2,dim> Y;



std::vector::active_cell_iterator> 
active_neighbors;

active_neighbors.reserve (GeometryInfo::faces_per_cell *

  GeometryInfo::max_children_per_face);



typename DoFHandler::active_cell_iterator 
cell_it(std_cxx11::get<0>(cell.iterators));





scratch_data.fe_midpoint_value.reinit (cell_it);



Tensor<1,dim> projected_gradient;



active_neighbors.clear ();



for (unsigned int face_no=0; 
face_no(cell.iterators)->at_boundary(face_no))

{



//just definition of abbreviation for the iterator to the 
face and the neighbor...

const typename DoFHandler::face_iterator

face = std_cxx11::get<0>(cell.iterators)->face(face_no);

const typename DoFHandler::cell_iterator

neighbor = std_cxx11::get<0
>(cell.iterators)->neighbor(face_no);





if (neighbor->active())

active_neighbors.push_back (neighbor);  //increase 
storage

else

{

// neighbor cell can be one level of coasered

// You should find all neighbor cell



// find a child of neighbor



if (dim == 1)  // special care is necessary for 1D case

{  



typename DoFHandler::cell_iterator 
neighbor_child = neighbor;

while (neighbor_child->has_children())

neighbor_child = neighbor_child->child 
(face_no==0 ? 1 : 0);



Assert (neighbor_child->neighbor(face_no==0 ? 1 : 0)

==std_cxx11::get<0
>(cell.iterators),ExcInternalError());



active_neighbors.push_back (neighbor_child);

}

else  // Most cases will fall under this case (2D and 
3D Problems..)

for (unsigned int subface_no=0; subface_no < 
face->n_children(); ++subface_no)

active_neighbors.push_back (std_cxx11::get<0
>(cell.iterators)->neighbor_child_on_subface(face_no,subface_no));

}

}



-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving 

[deal.II] Ideas on accuracy of Streamline diffusion method

[Jaekwang Kim]

Hi all, 

I would like ask today that possible factors that might affect, (badly) 
accuracy of streamline diffusion method for pure hyperbolic equation. 

I am solving hyperbolic equation , \nabla( beta \dot u ) = f(x), as it is 
suggested in tutorial step-9, 

Step-9 tutorial explained this very nicely, and I could get very 
well-behaved solution when I use linear element for my test and basis 
function. 

But, I am not receiving well behaved solution when I use higher order for 
my elements, the error in solution still has accuracy of O(h), not 
O(h^2)... 

I refered some references on Streamline Diffusion method, and papers tell 
me that I can expect at least p+(1/2) accuracy..but my accuracy is always 
O(h)... no matter the order of element I use. 

So, I suspect that I am being caught by somewhere that is only first order 
accurate. 

My weak form formation is exactly same with tutorial 9 



   // Advection Term

 scratch_data.fe_values.shape_value(i,q_point) * ( 
advection_directions[q_point] * 
scratch_data.fe_values.shape_grad(j,q_point)) //(beta grad u)



  // Stabilization term (Streamline Diffusion Method) 

+ *tau * * ( advection_directions[q_point] * 
scratch_data.fe_values.shape_grad(i,q_point) )

  *(advection_directions[q_point] * 
scratch_data.fe_values.shape_grad(j,q_point))





For tau , which determines the strength of diffusion term, was chosen to be 
'0.1h'. (Should this be chosen more wisely? Otherwise, I cannot expect 
'p+(1/2)' order for error? 


Another possible source of error might come from Integration Gauss 
quadrature, 

but the weak form (v,beta \cdot \nabla u) is locally '2p-1' polynomial, and 
I have used  enough Gauss Quadrature order for integrating the weak form.. 


Any one has idea on possible point I will have to check to attain (p+1/2) 
error behavior ? 



Always thank you, 


Regards, 


Jaekwang Kim 


 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: advection equation weak form in tutorial 9

[Jaekwang Kim]

Thanks, now I clearly understand what's going on!

Jaekwang 

On Monday, August 14, 2017 at 3:08:18 PM UTC-5, Jaekwang Kim wrote:
>
> Hi, all 
>
> I am little bit confused by weak form for advection problem written in 
> step-9 tutorial 
>
> in the its tutorial, in the middle of the line, it says that... 
>
> I am little bit unsure the red-marked form, 
>
> From the whole boundary integral, we put inflow case, to the right hand 
> side, and remaining part on left-hand-side is cases other than inflow 
> condition 
>
> so I thought, this should be.. 
>
>
> <https://lh3.googleusercontent.com/-_-hS916j4W8/WZICJB_jYWI/BEU/scEwE0R6BlYkOd4_vFYaAm-8rLuEOdqoACLcBGAs/s1600/Screen%2BShot%2B2017-08-14%2Bat%2B2.59.36%2BPM.png>
>  
>
> the left hand side term is boundary that is not inflow case... 
>
> I saw similar expression, in Dr. Bangerth video lecture 31, on hyperbolic 
> equation, here his weak form also marked 'not inflow boundary' remaining in 
> left hand side and inflow case, with boundary value description went to RHS 
> <https://lh3.googleusercontent.com/-cH5kYs776Is/WZICpHCSv8I/BEc/_xxauCWctgQjlxdx61-e9PCBuYZkLJINwCLcBGAs/s1600/Screen%2BShot%2B2017-08-14%2Bat%2B3.04.56%2BPM.png>
>
> . I am not sure i am just confusing some of indication... .. 
>
>
> thanks !
>
>
>
>
>
>
>
>
>
> This is a point we made in the introduction of step-3 
> <http://dealii.org/developer/doxygen/deal.II/step_3.html>. There, we 
> argued that to avoid this very kind of problem, one should get in the habit 
> of always multiplying with test functions *from the left* instead of from 
> the right to obtain the correct matrix right away. In order to obtain the 
> form of the linear system that we need, it is therefore best to rewrite the 
> weak formulation to 
>
> (vh+δβ⋅∇vh,β⋅∇uh)Ω*−(β⋅nvh,uh**)**∂**Ω**−*=(vh+δβ⋅∇vh,f)Ω−(β⋅nvh,g)∂Ω−
>
> and then to obtain
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] advection equation weak form in tutorial 9

[Jaekwang Kim]

Hi, all 

I am little bit confused by weak form for advection problem written in 
step-9 tutorial 

in the its tutorial, in the middle of the line, it says that... 

I am little bit unsure the red-marked form, 

>From the whole boundary integral, we put inflow case, to the right hand 
side, and remaining part on left-hand-side is cases other than inflow 
condition 

so I thought, this should be.. 


 

the left hand side term is boundary that is not inflow case... 

I saw similar expression, in Dr. Bangerth video lecture 31, on hyperbolic 
equation, here his weak form also marked 'not inflow boundary' remaining in 
left hand side and inflow case, with boundary value description went to RHS 


. I am not sure i am just confusing some of indication... .. 


thanks !









This is a point we made in the introduction of step-3 
. There, we argued 
that to avoid this very kind of problem, one should get in the habit of 
always multiplying with test functions *from the left* instead of from the 
right to obtain the correct matrix right away. In order to obtain the form 
of the linear system that we need, it is therefore best to rewrite the weak 
formulation to 

(vh+δβ⋅∇vh,β⋅∇uh)Ω*−(β⋅nvh,uh**)**∂**Ω**−*=(vh+δβ⋅∇vh,f)Ω−(β⋅nvh,g)∂Ω−

and then to obtain

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Can we say that the higher order method, the more accurate?

[Jaekwang Kim]

Dear Howe

I have solved this problem before, but I am not sure weather you are 
meeting same problem that I have encountered. 

If your mesh has unstructured structure, (not a rectangular one) you should 
higher order mapping when you asses finite element 
I remember that deal.ii default was Q1 mapping, so you might have been 
caught error form this... 

For example... I fixed my code as... 

Please look for 2 red lines I am impose same degree of mapping 

Thanks 

template 

void StokesProblem::initial_assemble_system ()

{

system_matrix=0;

system_rhs=0;

*const MappingQ mapping (degree);*

QGauss   quadrature_formula(degree+2);



FEValues fe_values (*mapping,* fe, quadrature_formula,

 update_values|

 update_quadrature_points  |

 update_JxW_values |

 update_gradients);



On Wednesday, August 9, 2017 at 3:06:39 AM UTC-5, Howe wrote:
>
> Dear Martin,
>
> Thanks for your rapid response.
> 1.
> The MappingQ is set to be the same as the order of velocity, as is shown 
> in the following code snippet:
>
>
>> *template *
>> *NS::NS (ParameterHandler )*
>> *:*
>> *   parameters (),*
>> *   degree (prm.get_integer("pressure degree")),*
>> *   fe( FE_Q(QGaussLobatto<1>(degree+2)), dim,*
>> *   FE_Q(QGaussLobatto<1>(degree+1)),   1),*
>> *   fe_scalar (FE_Q(QGaussLobatto<1>(degree+2))),*
>> *   dof_handler (triangulation),*
>> *   dof_handler_scalar (triangulation),*
>> *   mapping (degree+2),*
>> *   computing_timer (std::cout,*
>> *   TimerOutput::summary,**   TimerOutput::wall_times)*
>
>
> I am not quite sure whether the computation of the lift/drag in my code is 
> right, and my implementation  is almost the same as the one in this post: 
> https://groups.google.com/d/msg/dealii/rS6soTb69ig/C4QchAyEGwAJ 
> The only change is the first line :
>
> *QGauss face_quadrature_formula(degree+2);*
>
>
> 2. I am using deal.II 8.4.0 now, and I think i am not using manifold 
> description in my code
>
> 3. The mesh is shown as follows:
>
>
> <https://lh3.googleusercontent.com/-Vs0s6So28js/WYrCbyok9mI/BSc/rfs-1D1TA8c71IQZpOEyOoU3C1idNr1DQCLcBGAs/s1600/Image.png>
>
>
> Best,
>
> Howe
>
>
> 在 2017年8月9日星期三 UTC+8下午3:25:20，Martin Kronbichler写道：
>>
>> Dear Howe,
>>
>> How did you run your simulation? From your picture, it appears that a 
>> higher order method is worse at higher degrees than a lower order method, 
>> which does not match with my experience. If that were the case, nobody 
>> would use high orders. However, you need to bring many pieces in place to 
>> really get to the benefit of the high order method for somewhat more 
>> complicated examples such as the flow around a cylinder. Here is a list of 
>> things to look at:
>>
>> - Do you use a high-order polynomial mapping MappingQ of the same or 
>> higher degree as the interpolation space? Do you use this mapping in all 
>> routines that evaluate quantities, such as the usual assembly, the 
>> computation of the lift/drag, and so on?
>> - Do you use a manifold description that extends into the domain? (Look 
>> into TransfiniteInterpolationManifold.) Without, you will not get more than 
>> third order convergence.
>> - Do you have a good mesh around the area of interest? Flows around 
>> cylinders tend to be really really sensitive to the mesh quality around the 
>> cylinder.
>>
>> For the Navier-Stokes equations around the cylinder, if everything is 
>> done right one gets significantly improved results in terms of accuracy 
>> over the number of degrees of freedom up to degree (6,5) 
>> (velocity,pressure). Beyond that picture is less clear. At least with the 
>> meshes that we tried in our group it was not worth to go beyond. You can 
>> have a look a our results in section 5.4 and Figs. 9 and 10 of this 
>> preprint:
>> https://arxiv.org/pdf/1706.09252.pdf
>>
>> Best,
>> Martin
>>
>> On 09.08.2017 09:01, Howe wrote:
>>
>> Dear Jaekwang  
>>
>> Have you solved this problem? If yes, Could pls share your solution with 
>> us?
>> I am simulating a steady state flow over a cylinder, and the drag/lift 
>> coefficient shows an unexpected trend of change as i increase the 
>> discretization order and refine the mesh.
>>
>>
>> <https://lh3.googleusercontent.com/-Gz5932Zt6e0/WYqvbO3P-4I/AAAABSM/6hBHn1FO5W0P7X3SPfQ4iyRaldDYzOt3QCLcBGAs/s1600/Image.

[deal.II] Datapostprocessor first derivative based on what?

[Jaekwang Kim]

Hi, all 
I have asked specific implementation of Datapostprocessor Class a week ago. 

Adding to that.. I suddenly become curious how this class take 'the first 
derivative' form the solution vector. 

For example, if I solve flow problem with Taylor hood element as step-22 
tutorial program, and take the first derivatives of the acquired solution 
velocity using Datapostprocessor Class, 
I could easily get my target quantities (derivatives) for all same dof 
points. 

However, I come up with the fact that my finite elements is only 
C0-continuous...at the vertex where two different elements meet. 
(My taylor element is degree 2 for velocity and degree 1 for pressure. )
If then, on those vertexes, how the Datapostprocessor decides which element 
to refer for the first, or second derivatives on that dof? 

Thank you

Jaekwang Kim 




-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] how to use post_processor for FE_System

[Jaekwang Kim]

Thanks!! 
By the way, can I ask another? 

Actually, the quantity that I'm going to calculate is \dot(gamma), 
More explicitly...I attach one figure 

<https://lh3.googleusercontent.com/-Af70CldxGb4/WXC1hU2HssI/BDY/kTM9a2-uFmsJ1B1vNEpYoVbiDOUjGl6MgCLcBGAs/s1600/Screen%2BShot%2B2017-07-20%2Bat%2B8.50.15%2BAM.png>
So I modified post processor as, 

template 

void Magnitude::compute_derived_quantities_vector (

   const std::vector<Vector<
double> > & uh,

   const 
std::vector<std::vector<Tensor<1, dim> > >   & duh,

   const 
std::vector<std::vector<Tensor<2, dim> > >   & /*dduh*/,

   const 
std::vector<Point > & /*normals*/,

   const 
std::vector<Point > & evaluation_points,

   std::vector<Vector 
>   & computed_quantities) const

{

//const unsigned int n_quadrature_points = uh.size();



for (unsigned int i=0; i<computed_quantities.size(); i++)



{ Assert(computed_quantities[i].size() == 1,

   ExcDimensionMismatch (computed_quantities[i].size(), 1));

Assert(uh[i].size() == 3, ExcDimensionMismatch (uh[i].size(), 3));



Tensor<2,dim> grad_u;

for (unsigned int d=0; d<dim; ++d)

grad_u[d] = duh[i][d];



std::cout << "grad_u: " << grad_u[0][0] << " , " << grad_u[0][1] << 
std::end; 
// to check whether it actually gets value 



const SymmetricTensor<2,dim> shear_rate_tensor = symmetrize 
(grad_u);

computed_quantities[q](0) = shear_rate_tensor * shear_rate_tensor + 
std::pow( (uh[q](0) /evaluation_points[q](0)),2) );



}

}


but I got problem in two ways

1. the 'grad_u' vector I actually receiving is all zero vector I 
checked and it always gives me (dudx,dudy)=(0,0)
2. I need to use p[0] value, and I believe 'evaluation_points[q](0)' is 
what I need. 
but when I extract value I am receiving error as...

*grad_u: 0 , 0*

*make[3]: *** [CMakeFiles/run] Segmentation fault: 11*

*make[2]: *** [CMakeFiles/run.dir/all] Error 2*

*make[1]: *** [CMakeFiles/run.dir/rule] Error 2*

*make: *** [run] Error 2*


*Thanks...*

*Jaekwang *
 


On Wednesday, July 19, 2017 at 6:02:14 PM UTC-5, Wolfgang Bangerth wrote:
>
> On 07/19/2017 05:00 PM, Jaekwang Kim wrote: 
> > /See the entry in the Frequently Asked Questions of deal.II (linked to 
> > from http://www.dealii.org/) for a lot more information on what this 
> > error means and how to fix programs in which it happens./ 
>
> Have you looked here? :-) 
>
> (The answer is that the `Magnitude` object needs to live at least as 
> long as the `DataOut` object. You should declare it *before* the 
> `DataOut` object.) 
>
> Cheers 
>   W. 
>
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] how to use post_processor for FE_System

[Jaekwang Kim]

Thanks, !

I tried in a way you advised. so, now i tried. 



template 

void Magnitude::compute_derived_quantities_vector (

   const std::vector<Vector<
double> > ,

   const 
std::vector<std::vector<Tensor<1, dim> > >   & /*duh*/,

   const 
std::vector<std::vector<Tensor<2, dim> > >   & /*dduh*/,

   const 
std::vector<Point > & /*normals*/,

   const 
std::vector<Point > & /*evaluation_points*/,

   std::vector<Vector 
>   _quantities) const

{

for (unsigned int i=0; i<computed_quantities.size(); i++)



{ Assert(computed_quantities[i].size() == 1,

   ExcDimensionMismatch (computed_quantities[i].size(), 1));

*Assert(uh[i].size() == **3, ExcDimensionMismatch (uh[i].size(), 
3));*

computed_quantities[i](0) = std::sqrt(uh[i](0)*uh[i](0) + uh[i](1
)*uh[i](1));

}

}


and I use it as ...


template 

void

StokesProblem::output_results (const unsigned int 
refinement_cycle)  const

{   [...]   

DataOut data_out;

data_out.attach_dof_handler (dof_handler);

data_out.add_data_vector (solution, solution_names,

  DataOut::type_dof_data,

  data_component_interpretation);



*// added*

*Magnitude magnitude;*

*data_out.add_data_vector (solution, magnitude);  //here solution 
is block-vector solution *

*// added - end ;*



data_out.build_patches ();





std::ostringstream filenameeps;

filenameeps << "solution-"<< Utilities::int_to_string 
(refinement_cycle, 2)<< ".vtk";



std::ofstream output (filenameeps.str().c_str());

data_out.write_vtk (output);

[...]

}



but now I receive a error message as...where I don't have clue what's wrong 
now...

do you have any idea of it ? 


*An error occurred in line <114> of file 
 in function*

*virtual dealii::Subscriptor::~Subscriptor()*

*The violated condition was: *

*counter == 0*

*The name and call sequence of the exception was:*

*ExcInUse (counter, object_info->name(), infostring)*

*Additional Information: *

*Object of class 9MagnitudeILi2EE is still used by 1 other objects.*


*(Additional information: )*


*See the entry in the Frequently Asked Questions of deal.II (linked to from 
http://www.dealii.org/) for a lot more information on what this error means 
and how to fix programs in which it happens.*


On Wednesday, July 19, 2017 at 5:23:36 PM UTC-5, Wolfgang Bangerth wrote:
>
> On 07/19/2017 04:09 PM, Jaekwang Kim wrote: 
> > 
> > because my dof_handler has doc for pressure, size does not match... 
> > 
> > how can I avoid this problem ? 
>
> You just need to pass the *entire* solution vector, not just the first 
> block. Otherwise, DataOut cannot find out which elements of the vector 
> correspond to the degrees of freedom enumerated in your DoFHandler. 
>
> You will also have to change this: 
>
> Magnitude::compute_derived_quantities_vector ([...]) 
> { 
> [...] 
> Assert(uh[i].size() == 2, ExcDimensionMismatch (uh[i].size(), 2)); 
> [...] 
>
> The size will now be 3 (or probably dim+1), even though you are only 
> using the zeroth and first component of uh in computing the magnitude. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] how to use post_processor for FE_System

[Jaekwang Kim]

Hi, all 

I was practicing using post_processor, referring step-29 and step 31. 

If I want to use post_processor for vector-valued problem, where my 
solution is (u,v,p)

and I want to compute magnitude of velocity,  |u^2+v^2| ignoring pressure 
parts. 

i wrote as step 29; 

template 

class Magnitude : public DataPostprocessorScalar

{

public:

Magnitude ();



virtual

void

compute_derived_quantities_vector (const std::vector   
,

   const std::vector > > ,

   const std::vector > > ,

   const std::vector   
,

   const std::vector   
_points,

   std::vector 
_quantities) const;

  

};



template 

Magnitude::Magnitude ()

:

DataPostprocessorScalar ("ShearRate", update_values)

{}

// Remember that this function may only use data for which the

// respective update flag is specified by


template 

void

Magnitude::compute_derived_quantities_vector (

  const 
std::vector ,

  const 
std::vector > >   & /*duh*/,

  const 
std::vector > >   & /*dduh*/,

  const 
std::vector & /*normals*/,

  const 
std::vector & /*evaluation_points*/,

  
std::vector   _quantities

  ) const

{

for (unsigned int i=0; i of file 
 in 
function*

*void dealii::DataOut_DoFData, 2, 
2>::add_data_vector(const VectorType &, const 
DataPostprocessor &) [VectorType = 
dealii::Vector]*

*The violated condition was: *

*vec.size() == dofs->n_dofs()*

*The name and call sequence of the exception was:*

*Exceptions::DataOut::ExcInvalidVectorSize (vec.size(), dofs->n_dofs(), 
dofs->get_triangulation().n_active_cells())*

*Additional Information: *

*The vector has size 238 but the DoFHandler object says that there are 274 
degrees of freedom and there are 24 active cells. The size of your vector 
needs to be either equal to the number of degrees of freedom, or equal to 
the number of active cells.*




because my dof_handler has doc for pressure, size does not match... 

how can I avoid this problem ?  


-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google 

[deal.II] Re: weak form for advection-diffusion equation

[Jaekwang Kim]

Thanks for the response! 

I've read some articles that the governing equation (advection + diffusion) 
is unstable when it has large Peclet number.
(Even though it has diffusion term.)

What I am concerning is that it is how to imply SUPG stabilization method 
in deal.ii., because my future problem will include advection-diffusion 
term

In SUPG, people are considering residual term* ( b dot nabla - c nabla^2 u 
-f )*

.. and I thought this is very closely related to step-9 formation. 

But I am not still sure to assemble residual term in my weak form... 


<https://lh3.googleusercontent.com/-GaCm8gKtkbo/WUG4-V3gUHI/BCY/Htb_Gj81w-Y2Rc0pRogCsETEgNeBn7njACLcBGAs/s1600/Screen%2BShot%2B2017-06-14%2Bat%2B5.28.08%2BPM.png>


2017년 6월 14일 수요일 오후 3시 40분 50초 UTC-5, Bruno Turcksin 님의 말:
>
> Jaekwang,
>
> in step-9 the test function is a little bit strange because if you use 
> continuous finite element to discretize the advection equation, the scheme 
> is unstable. So you basically add a diffusion term to the equation to 
> stabilize the scheme. However, your equation already has a diffusion term 
> so you may not need a stabilization term.
>
> Best,
>
> Bruno
>
> On Wednesday, June 14, 2017 at 3:33:18 PM UTC-4, Jaekwang Kim wrote:
>>
>>
>> <https://lh3.googleusercontent.com/-qsDNe2Fge5w/WUGPdbVCfRI/BCI/egHPb83UKLIjnTGxmMfHiI4ZKGgXq4FhgCLcBGAs/s1600/Screen%2BShot%2B2017-06-14%2Bat%2B2.32.45%2BPM.png>
>>
>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: weak form for advection-diffusion equation

[Jaekwang Kim]

Adding up to this, I was trying as...


scratch_data.fe_values.shape_grad_grad (j,q_point)


but keep receiving errors 

*/Users/kimjaekwang/repos/advection_diffusion/step9/step-9.cc:487:53: **error: 
**no member named 'shape_grad_grad' in*

*  'dealii::FEValues<2, 2>'*

*   - scratch_data.fe_values.shape_grad_grad 
(j,q_point)*

point)



-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Weak form for advection - diffusion equation

[Jaekwang Kim]

Since
since my question include some mathematical expression, I attached figure 
for this...





















 copy_data.cell_matrix(i,j) += (

   

   (advection_directions[q_point] * 
scratch_data.fe_values.shape_grad(j,q_point)   *

   
(scratch_data.fe_values.shape_value(i,q_point) +



delta *

 (advection_directions[q_point] 
*

  
scratch_data.fe_values.shape_grad(i,q_point)))



//add diffusion term - should I 
consider stabilization term for this one too? Lets see.


+scratch_data.fe_values.shape_grad(i,q_point) 
*  scratch_data.fe_values.shape_grad(j,q_point) 

//stabilization term for 
diffusion

// (delta)(nabla^2 u , beta . 
nabla v)

+delta * ( - 
scratch_data.fe_values.shape_grad(j,q_point) * 
scratch_data.fe_values.shape_grad(j,q_point) )

   * 
(advection_directions[q_point] * 
scratch_data.fe_values.shape_grad(i,q_point) )



) *

   
scratch_data.fe_values.JxW(q_point));


but, I am not receiving solution I expected (it has manufactured 
solution...)

Is this approach wrong? or the direction is just correct, but I am making 
mistake somewhere in my code 


Thanks 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Distribute constraints when using Workstream class

[Jaekwang Kim]

I just understand what you meant ! 

Thanks !


2017년 4월 6일 목요일 오후 7시 24분 1초 UTC-5, Wolfgang Bangerth 님의 말:
>
> On 04/06/2017 02:55 PM, Jaekwang Kim wrote: 
> > template 
> > 
> > void 
> > 
> > Step4::copy_local_to_global (constAssemblyCopyData _data) 
> > 
> > { 
> > 
> > for(unsignedinti=0; i<copy_data.local_dof_indices.size(); ++i) 
> > 
> > { 
> > 
> > for(unsignedintj=0; j<copy_data.local_dof_indices.size(); ++j) 
> > 
> > system_matrix.add (copy_data.local_dof_indices[i], 
> > 
> >copy_data.local_dof_indices[j], 
> > 
> >copy_data.cell_matrix(i,j)); 
> > 
> > system_rhs(copy_data.local_dof_indices[i]) += 
> copy_data.cell_rhs(i); 
> > 
> > 
> > 
> > } 
> > 
> > 
> > 
> >  *constraints.distribute_local_to_global 
> > 
> (copy_data.cell_matrix,copy_data.cell_rhs,copy_data.local_dof_indices,system_matrix,system_rhs);*
>  
>
> > 
> > *// I am suspicious on this part...* 
> > 
> > } 
>
> You copy things twice into the global matrix and rhs, but only once taking 
> into accounts the constraints. Take a look at step-32, for example, to see 
> how 
> this function is supposed to look. 
>
> Best 
>   W. 
>
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: is there a way to refine mesh only in one direction

[Jaekwang Kim]

Here, I attached simple case that I encounter same error message. 

Most of code is just same at that of step-5, but I am using mesh-file that 
is made from other software. (I attach my mesh file too)



At constructor, there is a issue. 



template 

Step5::Step5 () :

  triangulation (Triangulation::maximum_smoothing),* //If I don't have 
this line, the code works well with anisotropy refinement.  *

*However, if I use Maximum_smoothing for my triangulation, I run into 
same error message.* 

  fe (1),

  dof_handler (triangulation)

{}




2017년 3월 31일 금요일 오후 6시 8분 16초 UTC-5, Wolfgang Bangerth 님의 말:
>
> On 03/31/2017 01:30 PM, Jaekwang Kim wrote: 
> > 
> > I found that what was the problem. 
> > 
> > I was using 'maximum' Meshsmoothing which does not allow 
> > 
> > So, It was not a bug of deal.ii but I was using wrong smoothing method 
> > that does not go along with anisotropy, 
>
> Well, it is still a bug that the library did not handle this in a more 
> gracefully way, i.e., that you did not get a polite error that suggests 
> what is happening here and why. 
>
> So I would still be interested in having a small testcase that shows the 
> problem and that would allow us to look into this some more... 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
##
#  CMake script for the step-5 tutorial program:
##

# Set the name of the project and target:
SET(TARGET "step-5")

# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
#FILE(GLOB_RECURSE TARGET_SRC  "source/*.cc")
#FILE(GLOB_RECURSE TARGET_INC  "include/*.h")
#SET(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC}) 
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
SET(TARGET_SRC
  ${TARGET}.cc
  )

# Usually, you will not need to modify anything beyond this point...

CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)

FIND_PACKAGE(deal.II 8.4 QUIET
  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
  )
IF(NOT ${deal.II_FOUND})
  MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to 
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()

DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()


mesh_stretch.inp
Description: chemical/gamess-input
/* -
 *
 * Copyright (C) 1999 - 2015 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Author: Wolfgang Bangerth, University of Heidelberg, 1999
 */



#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 

using namespace dealii;



template 
class Step5
{
public:
  Step5 ();
  void run ();

private:
  void setup_system ();
  void assemble_system ();
  void solve ();
  void output_results (const unsigned int cycle) const;
  void refine_grid (const unsigned int cycle);

  void set_anisotropic_flags ();

  Triangulation   triangulation;
  const MappingQ1 mapping;
  FE_Qfe;
  DoFHandler  dof_handler

Re: [deal.II] Re: is there a way to refine mesh only in one direction

[Jaekwang Kim]

Thanks Bruno ! and Dr. Wolfgang! 

I found that what was the problem.

I was using 'maximum' Meshsmoothing which does not allow 

So, It was not a bug of deal.ii but I was using wrong smoothing method that 
does not go along with anisotropy, 

I just changed my Meshsmoothing Method, that allow anisotropy refinement. 

For my problem, anisotropy refinement still works very well!

<https://lh3.googleusercontent.com/-6GN5r9R5tag/WN6uEamD0vI/BAk/IhFUj4xxd8s7eNxLUMN6XNTOIkDDr9MxQCLcB/s1600/Picture1.png>


Thank you agains all of you 

Jaekwang


2017년 3월 31일 금요일 오후 2시 22분 37초 UTC-5, Bruno Turcksin 님의 말:
>
> Jaekwang, 
>
> 2017-03-30 17:39 GMT-04:00 Jaekwang Kim <jaekw...@gmail.com >: 
>
> > //This one is problematic 
> > triangulation.prepare_coarsening_and_refinement (); 
> It could be a bug in the library, most people don't use anisotropic 
> refinement. This function is not always necessary (and it is not used 
> in step-30). So what happens if you remove it? For us to really help 
> you, we need a small code that triggers the error. 
>
> Best, 
>
> Bruno 
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Timo Heister and I on a podcast, talking about deal.II

[Jaekwang Kim]

Nice interview. 
I really enjoyed the cast. 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] is there a way to refine mesh only in one direction

[Jaekwang Kim]

Hi all,

I am running deal.ii for study non-newtonian fluids. 

I have a mesh which is attached. 

I am using AMR, but, when I use just conventional function, 
"cell->refine_flag()" for refinement, 

the flagged cell divide into 4. 

but as my problem only need high resolution in r-direction, 
I want to turn off theta directional refinement. 

So, in red square... I attached my desired refinement way. 

will that be way to accomplish this? , If needed, I think I can try merge 
some cell after refinement... 

if I need to describe question. please just let me know

Thanks. 


Bests, 
Jaekwang Kim 

<https://lh3.googleusercontent.com/-QyyWHbmWMC4/WNwDRpHKeEI/BAQ/3NlDrVNVqM0KsSdDDdu4XgePeI846Zj1gCLcB/s1600/Picture1.png>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: AMR , how to pass solution vector to refined mesh

[Jaekwang Kim]

Thank you!!
Just got what've meant. 
and fixed my code properly

Thanks ! 

Jaekwang Kim 

2017년 3월 21일 화요일 오후 12시 55분 7초 UTC-4, Jean-Paul Pelteret 님의 말:
>
> Hi Jaekwang,
>
> No, you've misunderstood in that I was just showing that to you as an 
> example. You don't need to use that exact Vector class. You were nearly at 
> the answer with what you've posted - you just need to read the compiler 
> error more carefully to understand what its telling you is wrong. The 
> deal.II BlockVector class is templatised on a number type, so must pass the 
> arguments along with the vector type when defining the vector template 
> argument of the SolutionTransfer class. So, in your case, this would be
>
> SolutionTransfer<dim,BlockVector > solution_transfer(dof_handler);
>>
>
> This I suspect will fix the problem you've shown above.
>
> Regards,
> Jean-Paul 
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: AMR , how to pass solution vector to refined mesh

[Jaekwang Kim]

thank you for your response. 
I am trying to following this... 

but it seems that I have to use the class, 
* SolutionTransfer<dim,TirilinoisWarppers::MPI::BlockVector>* to  transfer 
block vector 

Without using MPI, isn't it possible to assess this class? like... 
SolutionTransfer<dim,BlockVector> , while it shows error ,

*error: **use of*

*  class template 'BlockVector' requires template arguments*

*SolutionTransfer<dim,BlockVector> 
solution_transfer(dof_handler)...*

* ^~~*

*/Users/kimjaekwang/Programs/dealii/include/deal.II/n*

For now I don't have mpi in my computer.If there is no way to do this 
without using MPI, i have to install..


Thank you. 



2017년 3월 21일 화요일 오전 4시 25분 37초 UTC-4, Jean-Paul Pelteret 님의 말:
>
> Dear Jaekwang,
>
> Step-31 <https://www.dealii.org/8.4.1/doxygen/deal.II/step_31.html> 
> demonstrates 
> how to use the SolutionTransfer class with BlockVectors. In short, you need 
> to specify the vector type as a template argument to the class. Here's how 
> its done in that tutorial:
>
> SolutionTransfer<dim,TrilinosWrappers::MPI::Vector> 
> <https://www.dealii.org/8.4.1/doxygen/deal.II/classSolutionTransfer.html>
> temperature_trans(temperature_dof_handler); // non-block vector
> SolutionTransfer<dim,TrilinosWrappers::MPI::BlockVector> 
> <https://www.dealii.org/8.4.1/doxygen/deal.II/classSolutionTransfer.html>
> stokes_trans(stokes_dof_handler); // block vector
> triangulation.prepare_coarsening_and_refinement 
> <https://www.dealii.org/8.4.1/doxygen/deal.II/classTriangulation.html#ab9fa3177e0e43ab0cf243215d284a35a>
> ();
> temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);
> stokes_trans.prepare_for_coarsening_and_refinement(x_stokes);
>
> I hope this helps,
> Jean-Paul
>
> On Tuesday, March 21, 2017 at 12:26:59 AM UTC+1, Jaekwang Kim wrote:
>>
>> Dear, Dr. Bangerth
>>
>> Adding up to this, may I ask you one more question?
>>
>> After this,,, I was trying to transfer my block vector solution through 
>> Solutiontransfer class. but I am having difficulty in this. 
>>
>> For example, I have fluid problem vector value solution (u_x, u_y, 
>> pressure)
>>
>> at refie_mesh function, I tried...
>>
>> template 
>>
>> void
>>
>> StokesProblem::refine_mesh () //AMR based on VELOCITY GRADIENT
>>
>> {
>>
>>
>>  const FEValuesExtractors::Vector velocities (0);
>>
>> 
>>
>>  Vector estimated_error_per_cell 
>> (triangulation.n_active_cells());
>>
>> 
>>
>>  KellyErrorEstimator::estimate (dof_handler,QGauss> >(degree+1),typename FunctionMap::type(),
>>
>>  solution,estimated_error_per_cell, 
>> fe.component_mask(velocities));
>>
>>  
>>
>>  GridRefinement::refine_and_coarsen_fixed_number 
>> (triangulation,estimated_error_per_cell,0.3, 0.0);
>>
>> 
>>
>> triangulation.prepare_coarsening_and_refinement ();
>>
>>
>>
>> SolutionTransfer solution_transfer(dof_handler); *// Here 
>> my doc hander might include all (u,v,p) solution *
>>
>>
>>
>>
>>
>> After this I reinitialized my previous solution, where my current 
>> solution will be interpolated, for newly refined dof_handler. 
>>
>> And When I tried to transfer my solution 
>>
>>
>> 
>> solution_transfer.prepare_for_coarsening_and_refinement(solution.block(0
>> ));
>>
>> *//since this function only takes vector type inputs... I assume 
>> that I cannot transfer whole solution (u,v,p) at once? *
>>
>> 
>>
>> triangulation.execute_coarsening_and_refinement();
>>
>> 
>>
>> 
>>
>> //Transfer solution
>>
>> solution_transfer.interpolate(solution.block(0), 
>> previous_solution.block(0));
>>
>> 
>>
>> 
>> solution_transfer.prepare_for_coarsening_and_refinement(solution.block(1
>> ));
>>
>> 
>>
>> solution_transfer.interpolate(solution.block(1), 
>> previous_solution.block(1));
>>
>>
>> but I run into this error. 
>> I think it is because there two lines are conflicting  because they 
>> don't have same number of dof_handler...
>>
>> 1. SolutionTransfer solution_transfer(dof_handler);
>>
>> 2. 
>> solution_transfer.

Re: [deal.II] Re: AMR , how to pass solution vector to refined mesh

[Jaekwang Kim]

Dear. Dr, Bangerth

Adding up to this, I want to ask one more question. 

I was trying to transfer 'vector-valuled' solution (Block vector type) as I 
did for scalar solution (Vector type).

but it seems that the command line 

*solution_transfer.prepare_for_coarsening_and_refinement(solution); *


*error: **no matching member*

*  function for call to 'prepare_for_coarsening_and_refinement'*

*solution_transfer.prepare_for_coarsening_and_refinement(solution);*


is only able to transfer scalar solution. since I am receiving error 
message. 


I think I have to transfer my fluid problem solution (u,v,p) each by each 
to refined mesh. 


How can I do this...? (Is there a short way that I can assess each solution 
vector that live in BlockVector ?


Always thank you. 


Jaekwang Kim  

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: AMR , how to pass solution vector to refined mesh

[Jaekwang Kim]

Thank you ! 

You just pointed out exact problem that my code had!

I was not taking into constraints when I make new sparsity pattern! 


*  DoFTools::make_sparsity_pattern(dof_handler,*

*dsp,*

*constraints,*

*/*keep_constrained_dofs = */ false);*


2017년 3월 20일 월요일 오후 1시 3분 50초 UTC-4, Wolfgang Bangerth 님의 말:
>
> On 03/20/2017 10:58 AM, Jaekwang Kim wrote: 
> > 
> > I tried to assert assemble function as you suggested, but I didn't 
> > receive any error message. 
> > 
> > What make me more confuse me is the code is working with global 
> refinement. 
>
> Then are you building the sparsity pattern taking into account hanging 
> node constraints? Take a look how that is done in step-6, for example. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: AMR , how to pass solution vector to refined mesh

[Jaekwang Kim]

Thanks for your response. 

I tried to assert assemble function as you suggested, but I didn't receive 
any error message. 

What make me more confuse me is the code is working with global refinement.

if I input 100% cells to be refined.. 


 GridRefinement::refine_and_coarsen_fixed_number (triangulation,

 
estimated_error_per_cell,

 *1.0, 0.0);*


if I input 90 cells to be refined, it goes at least one cycle, 

but same or less than 80 cell, I keep receiving error message that... 

Might AMR be related to my program?  

 


*An error occurred in line <1668> of file 

 
in function*

*void dealii::internals::dealiiSparseMatrix::add_value(const LocalType, 
const size_type, const size_type, SparseMatrixIterator &) 
[SparseMatrixIterator = dealii::SparseMatrixIterators::Iterator, LocalType = double]*

*The violated condition was: *

*matrix_values->column() == column*

*The name and call sequence of the exception was:*

*typename SparseMatrix::ExcInvalidIndex(row, column)*

*Additional Information: *

*You are trying to access the matrix entry with index <1,305>, but this 
entry does not exist in the sparsity pattern of this matrix.*


*The most common cause for this problem is that you used a method to build 
the sparsity pattern that did not (completely) take into account all of the 
entries you will later try to write into. An example would be building a 
sparsity pattern that does not include the entries you will write into due 
to constraints on degrees of freedom such as hanging nodes or periodic 
boundary conditions. In such cases, building the sparsity pattern will 
succeed, but you will get errors such as the current one at one point or 
other when trying to write into the entries of the matrix.*





template 

void nonlinear::refine_grid ()

{

Vector estimated_error_per_cell (triangulation.n_active_cells());

KellyErrorEstimator::estimate (dof_handler,

QGauss(3),

typename FunctionMap::type(),

solution,

estimated_error_per_cell);

GridRefinement::refine_and_coarsen_fixed_number (triangulation,

 
estimated_error_per_cell,

 1.0, 0.0);



/* You need to tranfer solution (n-1 step) domain into previous 
solution (n step)*/

   

triangulation.prepare_coarsening_and_refinement ();

SolutionTransfer solution_transfer(dof_handler);

solution_transfer.prepare_for_coarsening_and_refinement(solution);

triangulation.execute_coarsening_and_refinement(); // I think n_dof() 
has increased here. You need to check this if necessary

dof_handler.distribute_dofs(fe);



   // Vector tmp(dof_handler.n_dofs()); //tmp - n step dof



previous_solution.reinit (dof_handler.n_dofs());



solution_transfer.interpolate(solution, previous_solution);  // 
interpolate (input, output)





DynamicSparsityPattern dsp(dof_handler.n_dofs());

DoFTools::make_sparsity_pattern (dof_handler, dsp);

sparsity_pattern.copy_from(dsp);



solution.reinit (dof_handler.n_dofs());

system_matrix.reinit (sparsity_pattern);

system_rhs.reinit (dof_handler.n_dofs());



constraints.clear ();

DoFTools::make_hanging_node_constraints (dof_handler,

 constraints);

VectorTools::interpolate_boundary_values (dof_handler,1,

  BoundaryValues(),

  constraints);

constraints.close ();



}

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: AMR , how to pass solution vector to refined mesh

[Jaekwang Kim]

Thanks for the response! 

I tried to put "Assert (system_matrix.n() == dof_handler.n_dofs(), 
ExcInternalError()); " when I assemble system , but i get the same error 
message, regardless of it 

*An error occurred in line <1668> of file 

 
in function*

*void dealii::internals::dealiiSparseMatrix::add_value(const LocalType, 
const size_type, const size_type, SparseMatrixIterator &) 
[SparseMatrixIterator = dealii::SparseMatrixIterators::Iterator, LocalType = double]*

*The violated condition was: *

*matrix_values->column() == column*

*The name and call sequence of the exception was:*

*typename SparseMatrix::ExcInvalidIndex(row, column)*

*Additional Information: *

*You are trying to access the matrix entry with index <3,13>, but this 
entry does not exist in the sparsity pattern of this matrix.*




By the way, i found that the code is working well when I refine my mesh 
globally, saying 100% of refinement. 



GridRefinement::refine_and_coarsen_fixed_number (triangulation,

 
estimated_error_per_cell,

 1.0, 0.0); 




it goes 5-refinement cycle perfectly. 


*Cycle 0:*

*   Number of active cells:   4*

*   Number of degrees of freedom: 25*

*dof_hander._n_dofs() : 25*

* Assemble system finished*

*   ||u_k-u_{k-1}|| = 9.91167e-06   Iteration Number: 52*

*   L2_error : 0.00233134*

*Cycle 1:*

*   Number of active cells:   16*

*dof_hander._n_dofs() : 81*

*Set up system finished*

* Assemble system finished*

*   ||u_k-u_{k-1}|| = 9.41577e-06   Iteration Number: 15*

*   L2_error : 0.000299119*

*Cycle 2:*

*   Number of active cells:   16*

*dof_hander._n_dofs() : 81*

*Set up system finished*

* Assemble system finished*

*   ||u_k-u_{k-1}|| = 9.98969e-07   Iteration Number: 2*

*   L2_error : 0.000287962*

*Cycle 3:*

*   Number of active cells:   16*

*dof_hander._n_dofs() : 81*

*Set up system finished*

* Assemble system finished*

*   ||u_k-u_{k-1}|| = 6.38877e-06   Iteration Number: 1*

*   L2_error : 0.000287717*

*[100%] Built target run*



However, if only 90% of cell is being refined, I could only go 3 cycle, 

and if then same of less than 80 I am receiving same errors... 

Might AMR be related to this...?



GridRefinement::refine_and_coarsen_fixed_number (triangulation,

   
  estimated_error_per_cell,

 0.8, 0.0); 



-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: AMR , how to pass solution vector to refined mesh

[Jaekwang Kim]

Thank you, Dr. Bangerth!

Step-15 tutorial was really helpful for me and I tried to used solution 
transfer class. 

I refined grid as follow. 

Because I am using Picard iteration, I use it as a previous_solution vector 
(in first refined mesh) when it is interpolated into refined meshes. 



template 

void nonlinear::refine_grid ()

{

Vector estimated_error_per_cell (triangulation.n_active_cells());

KellyErrorEstimator::estimate (dof_handler,

QGauss(3),

typename FunctionMap::type(),

solution,

estimated_error_per_cell);

GridRefinement::refine_and_coarsen_fixed_number (triangulation,

 
estimated_error_per_cell,

 0.3, 0.0);



/* You need to tranfer solution (n-1 step) domain into previous 
solution (n step)*/

   

triangulation.prepare_coarsening_and_refinement ();

SolutionTransfer solution_transfer(dof_handler);

solution_transfer.prepare_for_coarsening_and_refinement(solution);

triangulation.execute_coarsening_and_refinement(); // I think n_dof() 
has increased here. You need to check this if necessary

dof_handler.distribute_dofs(fe);



   // Vector tmp(dof_handler.n_dofs()); //tmp - n step dof



previous_solution.reinit (dof_handler.n_dofs());



solution_transfer.interpolate(solution, previous_solution);  // 
interpolate (input, output)



solution.reinit (dof_handler.n_dofs());



}


Also My setup_system function is is as follow 



template 

void nonlinear::setup_system (const unsigned int refinement_cycle)

{

if (refinement_cycle==0)

{ dof_handler.distribute_dofs (fe);



std::cout << "   Number of degrees of freedom: "

<< dof_handler.n_dofs()

<< std::endl;



solution.reinit (dof_handler.n_dofs());

previous_solution.reinit (dof_handler.n_dofs());

system_rhs.reinit (dof_handler.n_dofs());





for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)//

{//

previous_solution(i)=1;

}





DynamicSparsityPattern dsp(dof_handler.n_dofs());

DoFTools::make_sparsity_pattern (dof_handler, dsp);

sparsity_pattern.copy_from(dsp);



system_matrix.reinit (sparsity_pattern);







constraints.clear ();

DoFTools::make_hanging_node_constraints (dof_handler,

 constraints);



VectorTools::interpolate_boundary_values (dof_handler,1,

  BoundaryValues(),

  constraints);



constraints.close ();



}

else

{



//You don't need to make previous solution. Or solution vector , 
(it will be done on refinement step)

//What you should do?  set boundary condition again.

DynamicSparsityPattern dsp(dof_handler.n_dofs());

DoFTools::make_sparsity_pattern (dof_handler, dsp);

sparsity_pattern.copy_from(dsp);



system_matrix.reinit (sparsity_pattern);

system_rhs.reinit (dof_handler.n_dofs());





constraints.clear ();

DoFTools::make_hanging_node_constraints (dof_handler,

 constraints);

VectorTools::interpolate_boundary_values (dof_handler,1,

  BoundaryValues(),

  constraints);

constraints.close ();



std::cout << "Set up system finished" << std::endl;

}



}


I think this two function do everything I needed... but when I run my code, 
I run into error message when I assemble system in second time. 
what might be the problem.?

Always thank you!!

Jaekwang Kim 

*An error occurred in line <1668> of file 

 
in function*

*void dealii::internals::dealiiSparseMatrix::add_value(const LocalType, 
const size_type, const size_type, SparseMatrixIterator &) 
[SparseMatrixIterator = dealii::SparseMatrixIterators::Iterator<double, 
false>, LocalType = double]*

*The violated condition was: *

*matrix_values->column() == column*

*The name and call sequence of the exception was:*

*typename SparseMatrix::ExcInvalidIndex(row, column)*

*Additional Information: *

*You are trying to access the matrix entry with index <0,30>, but this 
entry does not exist in the sparsity pattern of this matrix.*


*The most common cause for this problem is that you used a method to build 
the sparsity pattern that did not (completely) take into account all of the 
entr

[deal.II] AMR , how to pass solution vector to refined mesh

[Jaekwang Kim]

Hi, all. 

I wonder if it is possible to pass solution vector to refined-meshed

I am solving non-linear problem with iterative method. 

I'd like to use coarse mesh solution as a initial solution to find accurate 
solution fast. 

However, my present code sets up system again on every refinement cycle. 
So, I lose all solution that I had in coarse mesh. 

I think it is really helpful for me if I knew a way to pass the previous 
solution to as a initial solution for next iterative method. 

I think it might be complicated... first vector size of solution might 
different when mesh is refined

is there anyone who has idea on this ?

Thanks, 

Jaekwang Kim 

template 

void nonlinear::run ()

{

unsigned int cycle_control=4;



for (unsigned int cycle=0; cycle<cycle_control; ++cycle)

{

std::cout << "Cycle " << cycle << ':' << std::endl;



if (cycle == 0)

{



GridGenerator::hyper_cube (triangulation, 0.0, 1.0);

triangulation.refine_global (1);




for (typename Triangulation::active_cell_iterator

 cell=triangulation.begin_active();

 cell!=triangulation.end(); ++cell)

{

for (unsigned int f=0; f<GeometryInfo::faces_per_cell; 
++f)

{

if (cell->face(f)->at_boundary() == true)

{

if ( std::abs(cell->face(f)->center()[0]-0.5)< 1e-12 
)

{

cell->face(f)->set_boundary_id(10);

}



// Boundary faces

static const double tol = 1e-12;

if ( std::abs(cell->face(f)->center()[0]-0.0)< tol )

{

cell->face(f)->set_boundary_id(1);

}

if ( std::abs(cell->face(f)->center()[0]-1.0)< tol )

{

cell->face(f)->set_boundary_id(1);

}

if ( std::abs(cell->face(f)->center()[1]-0.0)< tol )

{

cell->face(f)->set_boundary_id(1);

}

if ( std::abs(cell->face(f)->center()[1]-1.0)< tol )

{

cell->face(f)->set_boundary_id(1);

}

}

}

}





std::ofstream out ("grid-1.ucd");

GridOut grid_out;

grid_out.write_ucd (triangulation, out);





}



else



refine_grid ();



std::cout << "   Number of active cells:   "<< 
triangulation.n_active_cells() << std::endl;



setup_system ();

int iteration=0;

 

assemble_system ();

solve ();   



Vector difference(dof_handler.n_dofs());

Vector temp_vector(dof_handler.n_dofs());

Vector temp_difference(dof_handler.n_dofs());



difference=solution; // std::cout << " ||u_k-u_{k-1}||" << 
difference.l2_norm() << std::endl;

previous_solution = solution ;



double omega = 0.5; //relxation control number



int success_step=0 ;



do{



assemble_system();

solve();





difference = solution;

difference -= previous_solution; //This is temporary step



  

temp_vector=solution;

temp_vector.add(omega,difference);



temp_difference = temp_vector;

temp_difference -= solution; //Calculate different again



if (temp_difference.l2_norm()< difference.l2_norm())  
//temp_difference 
= temp_vector - solution , difference=solution-previous_solution

{



success_step+=1;



difference=temp_vector;

difference-=solution ;

solution=temp_vector;

}





previous_solution=solution ;



iteration+=1;



std::ofstream outerror2("iter_error.dat",std::ios::app);

outerror2 << degree << " " << triangulation.n_active_cells() << " 
" << iteration << " " << difference.l2_

[deal.II] Re: any suggestion how to run same code repeatedly with different boundary condition ?

[Jaekwang Kim]

Thank you !! solve the problem with the second method! 

I really appreciate this!

Jaekwang Kim 

2017년 3월 4일 토요일 오후 1시 37분 40초 UTC-6, Jean-Paul Pelteret 님의 말:
>
> Hi Jaekwang,
>
> The error message is pointing you to the fact that when you call run() a 
> second time you end up trying to rebuild the triangulation on top of the 
> one that you have the first time run() was called. There are a few ways 
> that you can work around this. Here are two suggestions:
>
> 1. Create a new instance of the laplace_problem each time you change the 
> boundary value
>
>   {
>
>bvalue = 1;
>
>Step4<2> laplace_problem_2d;
>
>laplace_problem_2d.run ();
>   }
>
>   {  
>
>bvalue = 2;
>   Step4<2> laplace_problem_2d;
>
>laplace_problem_2d.run ();
>
>  }
>
>
> 2. Separate the run() and make_grid() calls.
>
>   {
>
>Step4<2> laplace_problem_2d;
>   laplace_problem_2d.make_grid ();
>
>  
>bvalue = 1;
>laplace_problem_2d.run ();//Remove make_grid() from run()
>
>  
>
>bvalue = 2;
>
>laplace_problem_2d.run ();
>
>  
>
>  }
>
> You could of change your boundary value in other way, such as in a for 
> loop or using a MultipleParameterLoop 
> <https://www.dealii.org/8.4.1/doxygen/deal.II/classMultipleParameterLoop.html>
> .
>
> I hope this helps you.
>
> Regards,
> Jean-Paul
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] any suggestion how to run same code repeatedly with different boundary condition ?

[Jaekwang Kim]

Hi all, I want to run same problem for different boundary condition values 
that will come globally... for example, 
if i want to run step-4 code with  different boundary condition, 


#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 

#include 


#include 

#include 

#include 


#include 


using namespace dealii;


double bvalue ;


template 

class Step4

{

public:

  Step4 ();

  void run ();


private:

  void make_grid ();

  void setup_system();

  void assemble_system ();

  void solve ();

  void output_results () const;


  Triangulation   triangulation;

  FE_Qfe;

  DoFHandler  dof_handler;


  SparsityPattern  sparsity_pattern;

  SparseMatrix system_matrix;


  Vector   solution;

  Vector   system_rhs;

};




template 

class RightHandSide : public Function

{

public:

  RightHandSide () : Function() {}


  virtual double value (const Point   ,

const unsigned int  component = 0) const;

};




template 

class BoundaryValues : public Function

{

public:

  BoundaryValues () : Function() {}


  virtual double value (const Point   ,

const unsigned int  component = 0) const;

};





template 

double RightHandSide::value (const Point ,

  const unsigned int /*component*/) const

{

  double return_value = 0.0;

  for (unsigned int i=0; iget_dof_indices (local_dof_indices);

  for (unsigned int i=0; i boundary_values;

  VectorTools::interpolate_boundary_values (dof_handler,

0,

BoundaryValues(),

boundary_values);

  MatrixTools::apply_boundary_values (boundary_values,

  system_matrix,

  

Re: [deal.II] Assemble Righthand Side for vector-valued problem

[Jaekwang Kim]

Dr. Bangerth 


Thank you, 
I just fixed up what was wrong... 
in most cases, it is usually easy problems. 

The manufactured solution of velocities has not satisfies the continuity 
equation. 
(i.e. manufactured solution does not satisfy, div.u=0)...

Jaekwang Kim  

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Assemble Righthand Side for vector-valued problem

[Jaekwang Kim]

Thank you for your advice. 

I would make the problem simpler. Do you get the right solution if you 
> have a constant viscosity? If you iterate, do you get the solution after 
> one iteration? 
>

Yes, I checked this. 
After one iteration, I get the solution after one iteration when constant 
viscosity case.
 

>
> I will also note that your manufactured solution is symmetric, but your 
> computational solution is not. This provides you with a powerful way 
> because you don't actually need to look at the convergence, it's enough 
> to check that on a *coarse* mesh the solution is *symmetric*. It may be 
> that solution is already non-symmetric on 1 cell, or on a 2x2 mesh -- in 
> which case you already know that something is wrong. The question then 
> is why that is so -- is your rhs correct, for example? 
>

After modified some of mistakes, now I have following result for pressure 
field.  
With same manufactured pressure solution...


  




<2x2 mesh> 





<4x4 mesh> 







For velocity field, I still have at least qualitatively good result as I 
posted previous comment.
Tho it does not converge with mesh refinement when it is compared to 
manufactured solution in L2_norm base. 


 
 

> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Assemble Righthand Side for vector-valued problem

[Jaekwang Kim]

Thanks, I understand what you meant!!
I fixed it. 

Add up this.. can I ask your intuition for my another problem?

I am using deal-ii to analyze Non-newtonian fluids, which varies viscosity. 
For example, my viscosity eta is function of (u) field. I am solving this 
with iteration method.. (Picard) 
After coding, I am testing my code with method of manufacture. 

For squared domain, I implemented all boundary condition with my assumed 
solution and impose source term also to numerically generate my assumed 
solution 
however, the problem I got is also follow 


While my velocity field is qualitatively similar to manufactured solution, 
the pressure field is not following the exact form as seen below figure. 
It is very strange that..I can approach u-soltution form while my 
pressure does not


In addition , any mesh refinement would not make more accurate solution 
now. 

If you were me, what would you check more to find out what is problem?

<https://lh3.googleusercontent.com/-mOzFBuSfPGc/WJt7ZXct4rI/A_E/0UZ-rvkW6HMaGLnDA0uqM6cO51203vTbwCLcB/s1600/problem.png>





2017년 2월 8일 수요일 오후 1시 40분 46초 UTC-6, Wolfgang Bangerth 님의 말:
>
> On 02/08/2017 12:35 PM, Jaekwang Kim wrote: 
> >*/_local_rhs(i) += fe_values.shape_value(i,q) *_/* 
> > 
> > */_rhs_values[q] *_/* 
> > 
> > */_fe_values.JxW(q);_/* 
> > 
>
> rhs_values[q] is a Vector. You need to say which component of it 
> you want to multiply with. (Hint: The correct component is 
> fe.system_to_component_index(i).first). 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Assemble Righthand Side for vector-valued problem

[Jaekwang Kim]

gnored: could not match 'Point' against 'Vector'

operator * (const OtherNumberfactor,

*^*

*/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:2635:1:
 
note: *

  candidate template ignored: could not match 'SymmetricTensor<rank, 
dim,

  type-parameter-0-2>' against 'const double'

operator * (const SymmetricTensor<rank,dim,Number> ,

*^*

*/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:2655:1:
 
note: *

  candidate template ignored: could not match 'SymmetricTensor' against

  'Vector'

operator * (const Numberfactor,

*^*

*/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:2709:1:
 
note: *

  candidate template ignored: could not match 'SymmetricTensor<rank, 
dim,

  type-parameter-0-2>' against 'const double'

operator * (const SymmetricTensor<rank,dim,Number> ,

*^*

*/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:2736:1:
 
note: *

  candidate template ignored: could not match 'SymmetricTensor' against

  'Vector'

operator * (const Number factor,

*^*

*/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:2772:1:
 
note: *

  candidate template ignored: could not match 'SymmetricTensor<rank, 
dim,

  double>' against 'const double'

operator * (const SymmetricTensor<rank,dim> ,

*^*

*/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:2791:1:
 
note: *

  candidate template ignored: could not match 'SymmetricTensor' against

  'Vector'

operator * (const double factor,

*^*

*/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:3086:1:
 
note: *

  candidate template ignored: could not match 'SymmetricTensor<2, dim,

  type-parameter-0-1>' against 'const double'

operator * (const SymmetricTensor<2,dim,Number> ,

*^*

1 error generated.



I would appreciate if someone tells me what I am doing wrong...

Thank you !

Jaekwang Kim 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Assemble system with variational right hand side

[Jaekwang Kim]

Thanks for the reply, 

I just wanted to make sure that whether I implemented it right, 
but I still cannot sure whether I am doing correctly. Let me be with 
more specific questions. 

First, Yes I am considering 2-dimensional flow.
if I writes a code that looks as previous, 

template 

  void

  RightHandSide::vector_value (const Point ,

Vector   ) const

  {

values(0) =  pow(p[0],2);  

values(1) =  pow(p[1],2);

   * values(2) =  0.0;  // I assumed this always zero. because f is 2 
dimensional force, and we will not have third component of f*. 

  }

1) whether my interpretation of 'values(2) is correct' - marked with red 
line 

2) If we considered a vector value test function phi=(*v,*q)= (v_x, v_y , 
q), the right hand side term in weak form should look like this 

RHS = (*v, f* ) at domain o mega. 

since both *v* and* f* is vector that has two direction , (x and y) , the 
calculation should be conducted in both direction 
but in the following code, it seems that I am multiplying only the first 
component v and f... 
Do I understand this appropriately?

In this sense, I don't understand my code well... cuz it looks like...

*const unsigned int component_i =*

*   fe.system_to_component_index(i).first;   <- Here I only said the first *
*component...*

* local_rhs(i) += fe_values.shape_value(i,q) **

* rhs_values[q](component_i) **

* fe_values.JxW(q);*


Shouldn't this be some thing like this?  

*const unsigned int component_i =*

*   fe.system_to_component_index(i).first; *

*const unsigned int component_j = *

*  fe.system_to_component_index(j).second; *

* local_rhs(i) += (fe_values.shape_value(i,q) **
*rhs_values[q](component_i)) *

*(fe_values.shape_value(i,q) **
*rhs_values[q](component_j)) *
* fe_values.JxW(q);* 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Assemble system with variational right hand side

[Jaekwang Kim]

Hi. I have a question on how we can impose external force term, or right 
hand side, in Stokes flow. 

My governing equation is same with that of step-22, stokes flow 

-nabla^2 u + grad p = f 

-grad u = 0 

but I want to impose 'f', 

To do this, I made function that gives me f value on each point. 
For example, if I want function of f thats looks like f=(x^2,y^2); 
I would make function like this...

  template 

  class RightHandSide : public Function

  {

public:

  RightHandSide () : Function(dim+1) {}


  virtual void vector_value (const Point ,

Vector   ) const;


  };


  template 

  void

  RightHandSide::vector_value (const Point ,

Vector   ) const

  {

values(0) =  pow(p[0],2);

values(1) =  pow(p[1],2);

values(2) =  0.0;

  }

After that I may consider this 'f' term when I assemble system. 
but I don't see a way to how I should try modify following read lines of 
the code from step-22 

  for (; cell!=endc; ++cell)

  {

fe_values.reinit (cell);

local_matrix = 0;

local_rhs = 0;


right_hand_side.vector_value_list(fe_values.get_quadrature_points(),

  rhs_values);


for (unsigned int q=0; q

[deal.II] Re: Postprocessing on Vector Valued Problem

[Jaekwang Kim]

You were right

Oh I just got the result! 
I didn't know that dof_handler can automatically matching the sequence of 
the cell. 
Thank you for your time

Jaekwang Kim 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Vector Value Problem Solution initialization

[Jaekwang Kim]

Hi all. 

I have a question on the solution initialization for iterative methods, in 
vector valued problem. 

For example, I want to solve 2D flow, where my solution consist of two 
block, one is for velocities and the other is for pressure. 

First, I construct my initial solution as below 



*-Declaration  *


...BlockVector previous_solution;

...


*-Block build*


 {

  BlockDynamicSparsityPattern dsp (2,2);


  dsp.block(0,0).reinit (n_u, n_u);

  dsp.block(1,0).reinit (n_p, n_u);

  dsp.block(0,1).reinit (n_u, n_p);

  dsp.block(1,1).reinit (n_p, n_p);


  dsp.collect_sizes();


  DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false
);

  sparsity_pattern.copy_from (dsp);

}


*-reinitialization *


previous_solution.reinit (2);

  previous_solution.block(0).reinit (n_u);

  previous_solution.block(1).reinit (n_p);

  previous_solution.collect_sizes ();

 


After that , I wanted to set my x-componet of velocities by following 
command lines. 



*VectorTools::interpolate(dof_handler,initialize_solution(),previous_solution.block(0));*



However this function does not work because my dof_handler has more degree 
of freedom. (Degree of freedom of all velocities components + pressure 
components)

So I met error message which says

 

The violated condition was: 

vec.size() == dof.n_dofs()

The name and call sequence of the exception was:

ExcDimensionMismatch (vec.size(), dof.n_dofs())

Additional Information: 

Dimension 162 not equal to 187




It seems that I have to extract dofs that are only for velocity component 
to initialize my u_x velocities. 

Is there any way to do this? 

Or is there any other way to initialize my particular solution components 
only? 


& initialize_solution function looks like as follow.


template 

class initialize_solution : public Function

{

public:



initialize_solution () : Function() {}



virtual double value (const Point   ,

  const unsigned int component = 0) const;



};





template 

double

initialize_solution::value (const Point  ,

  const unsigned int component) const

{

Assert (component < this->n_components,

ExcIndexRange (component, 0, this->n_components));



double x=p[0];

double y=p[1];

double height = 1e-3; // unit of [m]

double return_value;



return_value=std::sin(y*Pi/height);



if (component == 0)  // impose only in x direction

return return_value;

else

return 0;



}




-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Very basic question on setting boundary id

[Jaekwang Kim]

Problem Solved! Thank you!! 

2016년 11월 13일 일요일 오후 5시 25분 20초 UTC-6, Daniel Arndt 님의 말:
>
> Jaekwang,
>
> It seems that you are missing
> constraints.distribute(solution);
> after solving the linear system. The constrained dofs are set to zero 
> during assembly and 
> this call is necessary to get the correct values for these in the solution 
> as well. You might want to have a look at
> https://www.dealii.org/8.4.1/doxygen/deal.II/group__constraints.html or 
> step-6.
>
> Best,
> Daniel
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Very basic question on setting boundary id

[Jaekwang Kim]

Thank you for being kindness! 

I followed your suggestions but I couldn't resolve the problem yet. 

Or did I make other mistake in another places?

>From a slight modification from STEP-4 tutorial , I have been trying solve 
Poisson Equation. 
Actually, I have manufactured solution to see my error behavior.
But I am testing whether I can really impose dirihlet boundary condition on 
the domain boundary. 

I would appreciate if you can read my code.cuz I couldn't resolve the 
problem, and having result like this...
although it is the boundary that I impose value of '5'...but it is set at 
center


 


2016년 11월 13일 일요일 오후 2시 43분 46초 UTC-6, Jean-Paul Pelteret 님의 말:
>
> Yes, well this is going to give strange results because you're setting the 
> boundary ID of internal faces (x=0.5 is the centreline of your geometry). 
> Before you set any boundary ID's, you should always first check to see that 
> you're actually on a boundary. The first few tutorials do go through this, 
> but I've provided you with a worked example to demonstrate the best 
> practise.
>
> Regards,
> J-P
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 1999 - 2016 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Author: Wolfgang Bangerth, University of Heidelberg, 1999
 */



#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 

#include 

#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 
#include 
#include 
#include 


#define Power(x,y)  ( std::pow(x,y))
#define Cos(x) ( std::cos(x) )
#define Sin(x) ( std::sin(x) )
#define Sqrt(x) ( std::sqrt(x) )
#define Pi 3.1415926535897932384626
#define E 2.71828182845904523


using namespace dealii;


template 
class Step4
{
public:
  Step4 (const unsigned int degree);
~Step4 ();
  void run ();

private:
  void make_grid ();
  void setup_system();
  void assemble_system ();
  void solve ();
  void refine_grid ();
  void error_evaluation();
  void output_results () const;

  const unsigned int degree;

  Triangulation   triangulation;
  FE_Qfe;
  DoFHandler  dof_handler;


  ConstraintMatrix constraints;

  SparsityPattern  sparsity_pattern;
  SparseMatrix system_matrix;

  double error = 99;

  Vector   solution;
  Vector   system_rhs;
};



/* -
 *
 * To test code with Manufactured solution 
 *  u=exp(-x^2-y^2);
 *
 * -*/

template 
class Solution : public Function
{
public:

Solution () : Function() {}

virtual double value (const Point   ,
  const unsigned int component = 0) const;

};


template 
double //generating exact solution : u=exp(-x^2-y^2);
Solution::value (const Point  ,
  const unsigned int component) const
{
Assert (component < this->n_components,
ExcIndexRange (component, 0, this->n_components));

double x=p[0];
double y=p[1];
double return_value;
double S=1; //Determine Strength of Singularity
double z=std::pow(x,2)+std::pow(y,2);
return_value=exp( (-1)*z/S );
return return_value;

}


/* -
 *
 * To test code with Manufactured solution
 * u=exp(-x^2-y^2);
 *
 * I am trying to impose value of '5' to see whether I can impose Dirihlet BC properly
 * -*/
template 
class BoundaryValues : public Function
{
public:
BoundaryValues () : Function() {}
virtual double value (const Point   ,
  const 

[deal.II] Re: Very basic question on setting boundary id

[Jaekwang Kim]

Thank you for quick reply!!! 

I tried your suggestions but I couldn't resolve the problem.

The more weird thing is if I commands as..


GridGenerator::hyper_cube (triangulation, 0, 1);

triangulation.refine_global ();


 for (typename Triangulation::active_cell_iterator

 cell=triangulation.begin_active();

 cell!=triangulation.end(); ++cell)

{

for (unsigned int f=0; fface(f)->center()[0]-0.5)< 1e-12 )


{


cell->face(f)->set_all_boundary_ids(1);



}





}

}






2016년 11월 13일 일요일 오후 1시 35분 3초 UTC-6, Jean-Paul Pelteret 님의 말:
>
> Hi Jaekwang,
>
> With what you've shown I'm going to guess that your problem might simply 
> be related to numerical / round-off error associated with floating point 
> calculations. 
>
> if ( cell->face(f)->center()[0]==1 ) // The LHS value is a computed 
>> double, and the RHS value is an integer which will be cast to a double for 
>> comparison
>
>
> Take a look at the mesh-generation part in the introduction to step-2 
> 
>  for 
> how one best implements these types of comparisons. They're particularly 
> relevant when dealing with the geometry, especially in more complex 
> scenarios. Admittedly the step-2 example is for a circular geometry but the 
> same concept applies: 
>
> if ( std::abs(cell->face(f)->center()[0]) < 1e-10 ) ...
>
> if ( std::abs(cell->face(f)->center()[0] - 1.0) < 1e-10 ) ...
>
> You can almost never be certain that any calculated floating point value 
> will be an exact match to an implicitly cast integer value right down to 
> machine precision, which is what will be checked with the equality operator.
>>
>>
> I hope that this fixes your issue.
>
> Regards,
> J-P
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Very basic question on setting boundary id

[Jaekwang Kim]

Hi all, 

I am learning grid generator how to set boundary condition on given domain.
It is very basic, but I couldn't even set boundary id for squared mesh 

I want to generate squared computational domain and I want to give boundary 
id of 1 to all the outer boundary of the square. 

So I used following lines 




GridGenerator::hyper_cube (triangulation, 0, 1);

triangulation.refine_global ();



std::cout << "Number of active cells: "<< 
triangulation.n_active_cells()<< std::endl;







for (typename Triangulation::active_cell_iterator

 cell=triangulation.begin_active();

 cell!=triangulation.end(); ++cell)

{

for (unsigned int f=0; f<GeometryInfo::faces_per_cell; 
++f)

{





if ( cell->face(f)->center()[0]==0 )


{

cell->face(f)->set_all_boundary_ids(1);



}



if ( cell->face(f)->center()[0]==1 )



{

cell->face(f)->set_all_boundary_ids(1);



}



if ( cell->face(f)->center()[1]==0 )



{

cell->face(f)->set_all_boundary_ids(1);



}

if ( cell->face(f)->center()[1]==1 )



{

cell->face(f)->set_all_boundary_ids(1);



}





}


However, when I checked my boundary by giving arbitrary boundary value of 
 '5', the solution shape looks as follow which is obviously wrong..


Are my code lines wrong? I thought that the faces of cell is each edge of a 
cell, and examined whether their centers are x=0,1 or y=0,1. 

and If then, I set id on 1...



<https://lh3.googleusercontent.com/-toZV8W-J-lw/WCi4Kldx3KI/A88/OF-tsoWgJpILwNta9WdYtUEyxMx77EBsQCLcB/s1600/Screen%2BShot%2B2016-11-13%2Bat%2B12.59.38%2BPM.png>



Always thank you all you for fast response!


Jaekwang Kim 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: hp convergence for non-linear solver.

[Jaekwang Kim]

Thanks for replying me ! 

I moved to squared domain to avoid other possible error when I imply 
boundary condition, but I still have same manufactured solution e^(-x^2-y^2)
I have to say that I have some changed error curve after moving to new mesh 
domain. 
Which is slightly better, but still unsatisfying results. 

<https://lh3.googleusercontent.com/-UR8qKv0nh0o/WBuw-Lg7gFI/A8k/FznfGdlllE8OifxsvLJ_Omudg90D611EgCLcB/s1600/error_squared_domain.jpg>

To see error per cell, I referred step 27, and make estimated error data 
out. 

I used following code lines



 Vector estimated_error_per_cell (triangulation.n_active_cells());



KellyErrorEstimator::estimate (dof_handler,

QGauss(3),

typename FunctionMap::type(),

solution,

estimated_error_per_cell);


In fact, I don't know how kelly error estimator works...
but just checked that it usually works well with Poisson Equations.

<https://lh3.googleusercontent.com/-JVGEEcGSQ6E/WBuwfOpt-xI/A8g/b8Wjo_cqspAgpz0HEViKz7y6jCz_5ggJQCLcB/s1600/Solution%2Band%2BError%2BEstimator.jpg>

Seeing pictures above,,, when p=2 solution looks seemingly good, but 
estimated error_per_cell should be wrong. 
Since, exact solution the closer points to x=(0,0) (or left down corner in 
pic) has  sharp gradient, thus the larger error is expected. 
Think estimated error case for p=1 is correct. 

Now, I confirmed that the code is not working at higher degree. (p>=2)...
but it is not because of boundary condition...since I have four boundary 
(top,bottom,right,left...) but the points which show highest error is some 
what irrealted to boundary...

Do you have any idea about this? 

Always thank you for your time and advice!

Best regards,
Jaekwang Kim 



2016년 11월 2일 수요일 오후 8시 30분 19초 UTC-5, Wolfgang Bangerth 님의 말:
>
> On 11/02/2016 07:17 PM, Jaekwang Kim wrote: 
> > 
> > I'd like check the visualized error per cell. 
> > Is there any module in deal.ii that enables this...? 
> > 
>
> Yes, you just need to create a vector with as many entries as there are 
> cell 
> (e.g., the output of VectorTools::integrate_difference) and attach that to 
> the 
> DataOut object you have. step-27 shows how to do that, but it in essence 
> works 
> just like for vectors that have as many entries as there are DoFs. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: hp convergence for non-linear solver.

[Jaekwang Kim]

Thanks for replying! 

Could you explain more?

I considered mapping degree only when I assemble system... for example. 

template 

void msurface::assemble_system ()

{

  const MappingQ mapping (degree);

  const QGauss  quadrature_formula(degree+2);



  FEValues fe_values (mapping, fe, quadrature_formula,

   update_values|  update_gradients |

   update_quadrature_points  |  update_JxW_values);






but can I also input a non-default Mapping degree when I use 
Vectortools::integrate_difference function?

I tried as ... but it seems that its not a proper way to use this 
function...?

template 

void msurface::evaluate_error()

{

Vector difference_per_cell (triangulation.n_active_cells());



const MappingQ mapping (degree);



VectorTools::integrate_difference (degree, 
dof_handler,solution,Solution(),

   
difference_per_cell,QGauss(degree+2),VectorTools::L2_norm);


How can I use non-default Mapping in here? (I wonder specific lines of 
code...)

Thank for giving me a help!


Regards, 

Jaekwang Kim

2016년 11월 1일 화요일 오후 5시 19분 35초 UTC-5, Daniel Arndt 님의 말:
>
> Jaekwang, 
>
> [...]
>>
>> VectorTools::integrate_difference 
>> (dof_handler,solution,Solution(),
>>
>>difference_per_cell
>> *,QGauss(degree+**2**)*,VectorTools::L2_norm);
>>
> somegthing that immediately comes to mind is that you don't use a 
> non-default Mapping in here although your boundary is curved.
> Make sure that you use an appropriate Mapping everywhere!
>
> Best,
> Daniel
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: hp convergence for non-linear solver.

[Jaekwang Kim]

<https://lh3.googleusercontent.com/-JiteM_JBmDk/WBkSG1wvCXI/A78/dVw6vhrCvcYrgwozJcO-gy76rNOuXQYBgCLcB/s1600/case1_p2_80.jpg>
Thanks for replying! 

Here I attached solution for p=2, and it seems that it has right share. 

When you said 'right quadrature' dose that meanquadrature for 
integration when we assemble system? 

To attach some of my code. 

template 

void msurface::assemble_system ()

{

  const MappingQ mapping (degree);

  const QGauss  quadrature_formula(degree+2);



  FEValues fe_values (mapping, fe, quadrature_formula,

   update_values|  update_gradients |

   update_quadrature_points  |  update_JxW_values);





I am pretty sure that I am using right quadrature here 


but I am not sure whether I am using right quadrature when I compare my 
numerical solution to exact solution.. 


VectorTools::integrate_difference (dof_handler,solution,Solution(),

   difference_per_cell
*,QGauss(degree+**2**)*,VectorTools::L2_norm);



Thanks again for fast reply!! 




2016년 11월 1일 화요일 오후 5시 4분 40초 UTC-5, Bruno Turcksin 님의 말:
>
> Jaekwang,
>
> On Tuesday, November 1, 2016 at 5:43:36 PM UTC-4, Jaekwang Kim wrote:
>>
>> I think that I didn't make any mistake on my code lines because at least 
>> my numerical solution is converging for the case p=1. 
>>
>> but I couldn't get such convergence for the case p=2, and 3. I have tried 
>> to change tolerance number in iterative scheme, but it didn't work. 
>>
>>
>> Are there other things that I failed to consider? 
>>
>> Or did I make wrong way to test the code? 
>>
> What you did looks good. You should check that you use the right 
> quadrature everywhere. If somewhere in your code you use a low order 
> quadrature, it will work for p=1 but you will have a larger error when you 
> increase the quadrature. You should also look at the solution for p=2 and 
> p=3. This might help you understand what is going on.
>
> Best,
>
> Bruno
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] hp convergence for non-linear solver.

[Jaekwang Kim]

Hi, all. 

I have asked whether numerical errors also behave ~O(h^p) , even though 
governing equation is  non-linear and thus, a solver includes iterative 
method. 
I recall that Dr. Bangerth replied me it is "yes" if we imply small enough 
tolerance on iteration tolerance. 
I am trying to confirm this testing the code with Picard Iteration, ( Dr. 
Bangerth's video lecture from 31.5~31.7)

To bring the governing equation again 

<https://lh3.googleusercontent.com/-TYzXkQCHcTA/WAksZOt6tjI/A7E/kX0e1wXWcfMcKc-dDYeq2rOOmXi-YtHjwCLcB/s1600/Screen%2BShot%2B2016-10-20%2Bat%2B3.42.54%2BPM.png>

Since I want to test the code really has correct h-p convergence, I used 
the method of manufactured solution. 


To be specific, I assumed solution as 


" u(x,y)= 
<https://lh3.googleusercontent.com/-dBW9xWRIpWE/WBkKfMK0QPI/A7c/SURwfidSLLAUL0ZrA94wn_sheykuZTDyACLcB/s1600/Screen%2BShot%2B2016-11-01%2Bat%2B4.34.36%2BPM.png>
"


and this gives me a right hand side function, f ,


f= 
<https://lh3.googleusercontent.com/-YdYt94kXCtc/WBkK3hcHHTI/A7g/iLbRiSE6hJk9Vvdpkc9uhQFhumJOuzpCQCLcB/s1600/Screen%2BShot%2B2016-11-01%2Bat%2B4.36.10%2BPM.png>where
 
z=x^2+y^2



After that I imposed boundary condition as exact with u(x,y) = 
<https://lh3.googleusercontent.com/-dBW9xWRIpWE/WBkKfMK0QPI/A7c/SURwfidSLLAUL0ZrA94wn_sheykuZTDyACLcB/s1600/Screen%2BShot%2B2016-11-01%2Bat%2B4.34.36%2BPM.png>
 =g 


Lastly, I compared my numerical solution with manufactured exact solution 
using Vectortools::integrate_difference. l2norm


template 

void msurface::evaluate_error()

{

Vector difference_per_cell (triangulation.n_active_cells());



VectorTools::integrate_difference (dof_handler,solution,Solution(),

   
difference_per_cell,QGauss(degree+2),VectorTools::L2_norm);



const double L2_error = difference_per_cell.l2_norm();


std::cout << "   L2_error : " << L2_error << std::endl;

   

error=L2_error;

}


However, I the error graph I get was as follow...

<https://lh3.googleusercontent.com/-YGzF49KtLZ8/WBkL1TaMVfI/A7s/yYA3RR2pPP0gTEbjRPkqaviDfjZJXDP7gCLcB/s1600/case1.jpg>



I think that I didn't make any mistake on my code lines because at least my 
numerical solution is converging for the case p=1. 

but I couldn't get such convergence for the case p=2, and 3. I have tried 
to change tolerance number in iterative scheme, but it didn't work. 


Are there other things that I failed to consider? 

Or did I make wrong way to test the code? 



Always thank you for all smart guys here!! 

Jaekwang Kim  

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: How to use sacado package

[Jaekwang Kim]

Always thank you for the comments!!

I'd like to follow the first choice, installing MPI on my machine. 

I installed openmpi-2.0.1. 

Coming to the original problem, I tried install trillions 12.4.2 again. 
but now I received this error message while i am installing trilinos. 

what might be the problem now? .. Do I use other packages for mph? 

"
.

[  5%] Building CXX object 
packages/teuchos/core/src/CMakeFiles/teuchoscore.dir/Teuchos_VerbosityLevel.cpp.o

[  5%] Building CXX object 
packages/teuchos/core/src/CMakeFiles/teuchoscore.dir/Teuchos_Workspace.cpp.o

[  5%] Building CXX object 
packages/teuchos/core/src/CMakeFiles/teuchoscore.dir/Teuchos_dyn_cast.cpp.o

[  5%] Building CXX object 
packages/teuchos/core/src/CMakeFiles/teuchoscore.dir/Teuchos_stacktrace.cpp.o

[  5%] *Linking CXX shared library libteuchoscore.dylib*

Undefined symbols for architecture x86_64:

  "_MPI_Allgather", referenced from:

  Teuchos::GlobalMPISession::allGather(int, Teuchos::ArrayView 
const&) in Teuchos_GlobalMPISession.cpp.o

  "_MPI_Allreduce", referenced from:

  Teuchos::GlobalMPISession::sum(int) in Teuchos_GlobalMPISession.cpp.o

  "_MPI_Barrier", referenced from:

  Teuchos::GlobalMPISession::barrier() in Teuchos_GlobalMPISession.cpp.o

  "_MPI_Comm_rank", referenced from:

  Teuchos::GlobalMPISession::initialize(std::__1::basic_ostream*) in Teuchos_GlobalMPISession.cpp.o

  "_MPI_Comm_size", referenced from:

  Teuchos::GlobalMPISession::initialize(std::__1::basic_ostream*) in Teuchos_GlobalMPISession.cpp.o

  "_MPI_Finalize", referenced from:

  Teuchos::GlobalMPISession::~GlobalMPISession() in 
Teuchos_GlobalMPISession.cpp.o

  "_MPI_Get_processor_name", referenced from:

  Teuchos::GlobalMPISession::GlobalMPISession(int*, char***, 
std::__1::basic_ostream*) in 
Teuchos_GlobalMPISession.cpp.o

  "_MPI_Init", referenced from:

  Teuchos::GlobalMPISession::GlobalMPISession(int*, char***, 
std::__1::basic_ostream*) in 
Teuchos_GlobalMPISession.cpp.o

  "_MPI_Initialized", referenced from:

  Teuchos::GlobalMPISession::GlobalMPISession(int*, char***, 
std::__1::basic_ostream*) in 
Teuchos_GlobalMPISession.cpp.o

  Teuchos::GlobalMPISession::initialize(std::__1::basic_ostream*) in Teuchos_GlobalMPISession.cpp.o

  Teuchos::Time::Time(std::__1::basic_string const&, bool) in 
Teuchos_Time.cpp.o

  Teuchos::Time::start(bool) in Teuchos_Time.cpp.o

  Teuchos::Time::Time(std::__1::basic_string const&, bool) in 
Teuchos_Time.cpp.o

  Teuchos::Time::wallTime() in Teuchos_Time.cpp.o

  Teuchos::Time::stop() in Teuchos_Time.cpp.o

  ...

  "_MPI_Wtime", referenced from:

  Teuchos::Time::Time(std::__1::basic_string const&, bool) in 
Teuchos_Time.cpp.o

  Teuchos::Time::start(bool) in Teuchos_Time.cpp.o

  Teuchos::Time::Time(std::__1::basic_string const&, bool) in 
Teuchos_Time.cpp.o

  Teuchos::Time::wallTime() in Teuchos_Time.cpp.o

  Teuchos::Time::stop() in Teuchos_Time.cpp.o

  Teuchos::Time::totalElapsedTime(bool) const in Teuchos_Time.cpp.o

  "_ompi_mpi_comm_world", referenced from:

  Teuchos::GlobalMPISession::initialize(std::__1::basic_ostream*) in Teuchos_GlobalMPISession.cpp.o

  Teuchos::GlobalMPISession::barrier() in Teuchos_GlobalMPISession.cpp.o

  Teuchos::GlobalMPISession::sum(int) in Teuchos_GlobalMPISession.cpp.o

  Teuchos::GlobalMPISession::allGather(int, Teuchos::ArrayView 
const&) in Teuchos_GlobalMPISession.cpp.o

  "_ompi_mpi_int", referenced from:

  Teuchos::GlobalMPISession::sum(int) in Teuchos_GlobalMPISession.cpp.o

  Teuchos::GlobalMPISession::allGather(int, Teuchos::ArrayView 
const&) in Teuchos_GlobalMPISession.cpp.o

  "_ompi_mpi_op_sum", referenced from:

  Teuchos::GlobalMPISession::sum(int) in Teuchos_GlobalMPISession.cpp.o

ld: symbol(s) not found for architecture x86_64

clang: error: linker command failed with exit code 1 (use -v to see 
invocation)

make[2]: *** [packages/teuchos/core/src/libteuchoscore.12.4.2.dylib] Error 1

make[1]: *** [packages/teuchos/core/src/CMakeFiles/teuchoscore.dir/all] 
Error 2

make: *** [all] Error 2






2016년 10월 22일 토요일 오후 9시 14분 1초 UTC-5, Wolfgang Bangerth 님의 말:
>
>
> Jaekwang, 
>
> > 
> */Users/jaekwangjk/Programs/trilinos-12.4.2-Source/packages/teuchos/core/src/Teuchos_GlobalMPISession.cpp:49:12:
>  
>
> > **fatal error: * 
> > 
> > *  'mpi.h' file not found* 
>
> If you configure Trilinos to use MPI (which you did when you called 
> Trilinos' 
> cmake script), then you obviously have to have MPI installed 

Re: [deal.II] Re: How to use sacado package

[Jaekwang Kim]

Thanks for the comment. 

I'd like to follow the first option to install MPI on my system. 

It seems that there are different options that I can choose. 

Ready page in deal.ii linked me to the site "http://mpi-forum.org; and it 
seems that there are different packages of mpi. 
Does it matter to choose any of packages? 

For example, I am going to install Open MPI newest version of  2.0.1?

>
> Jaekwang, 
>
> > 
> */Users/jaekwangjk/Programs/trilinos-12.4.2-Source/packages/teuchos/core/src/Teuchos_GlobalMPISession.cpp:49:12:
>  
>
> > **fatal error: * 
> > 
> > *  'mpi.h' file not found* 
>
> If you configure Trilinos to use MPI (which you did when you called 
> Trilinos' 
> cmake script), then you obviously have to have MPI installed on your 
> system. 
> You have two options: You can either install the MPI package on your 
> machine, 
> or you can tell Trilinos not to use MPI. That choice is yours. 
>
> Best 
>   W. 
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to use sacado package

[Jaekwang Kim]

Thanks Daniel! 

I tried to make install Trillions , but after the command "make install" I 
run into other error. 

*/Users/jaekwangjk/Programs/trilinos-12.4.2-Source/packages/teuchos/core/src/Teuchos_GlobalMPISession.cpp:49:12:
 
**fatal error: *

*  'mpi.h' file not found*

#  include "mpi.h"


. 

Is this because I am not uniting cluster?

If then,, Cant I uses Trilinos package?


Thank you for your kindness!

2016년 10월 22일 토요일 오전 4시 2분 50초 UTC-5, Daniel Arndt 님의 말:
>
> Jaekwang,
>
> It seems that no fortran compiler is found. You should be able to 
> configure Trilinos without fortran using
> -DTrilinos_ENABLE_Fortran:BOOL=OFF
>
> After that it's the usual make install [1].
>
> Best,
> Daniel
>
> [1] https://www.dealii.org/developer/external-libs/trilinos.html
>
> Am Samstag, 22. Oktober 2016 01:44:26 UTC+2 schrieb Jaekwang Kim:
>>
>>
>> I'd like to use sacado package for automatic differentiation. 
>>
>>
>> However, I have installed deal.ii on my computer with trilinos 
>> configuration set off. 
>>
>>
>> I looked back deal.ii installation manual and got to know that I have to 
>> install trilnos package first. 
>>
>>
>> However, when I try to install trilnos following instruction, 
>>
>>
>> cd trilinos-12.4.2
>> mkdir build
>> cd build
>>
>> cmake\
>> -DTrilinos_ENABLE_Amesos=ON  \
>> -DTrilinos_ENABLE_Epetra=ON  \
>> -DTrilinos_ENABLE_Ifpack=ON  \
>> -DTrilinos_ENABLE_AztecOO=ON \
>> -DTrilinos_ENABLE_Sacado=ON  \
>> -DTrilinos_ENABLE_Teuchos=ON \
>> -DTrilinos_ENABLE_MueLu=ON   \
>> -DTrilinos_ENABLE_ML=ON  \
>> -DTrilinos_VERBOSE_CONFIGURE=OFF \
>> -DTPL_ENABLE_MPI=ON  \
>> -DBUILD_SHARED_LIBS=ON   \
>> -DCMAKE_VERBOSE_MAKEFILE=OFF \
>> -DCMAKE_BUILD_TYPE=RELEASE   \
>> -DCMAKE_INSTALL_PREFIX:PATH=$HOME/share/trilinos \
>> ../
>>
>> make install
>>
>>
>> I got error message as follows...
>>
>>
>>
>>
>> "
>>
>> The Fortran compiler identification is unknown
>>
>> CMake Error at 
>> cmake/tribits/core/package_arch/TribitsGlobalMacros.cmake:1638 
>> (ENABLE_LANGUAGE):
>>
>>   No CMAKE_Fortran_COMPILER could be found.
>>
>>
>>   Tell CMake where to find the compiler by setting either the environment
>>
>>   variable "FC" or the CMake cache entry CMAKE_Fortran_COMPILER to the 
>> full
>>
>>   path to the compiler, or to the compiler name if it is in the PATH.
>>
>> Call Stack (most recent call first):
>>
>>   cmake/tribits/core/package_arch/TribitsProjectImpl.cmake:188 
>> (TRIBITS_SETUP_ENV)
>>
>>   cmake/tribits/core/package_arch/TribitsProject.cmake:93 
>> (TRIBITS_PROJECT_IMPL)
>>
>>   CMakeLists.txt:89 (TRIBITS_PROJECT)
>>
>> "
>>
>>
>> What's the problem? 
>>
>> and I would very thank if someone tell me how to change deal.ii 
>> configuration to set on Trilinos after that
>>
>>
>> Thank you! 
>>
>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] How to use sacado package

[Jaekwang Kim]

I'd like to use sacado package for automatic differentiation. 


However, I have installed deal.ii on my computer with trilinos 
configuration set off. 


I looked back deal.ii installation manual and got to know that I have to 
install trilnos package first. 


However, when I try to install trilnos following instruction, 


cd trilinos-12.4.2
mkdir build
cd build

cmake\
-DTrilinos_ENABLE_Amesos=ON  \
-DTrilinos_ENABLE_Epetra=ON  \
-DTrilinos_ENABLE_Ifpack=ON  \
-DTrilinos_ENABLE_AztecOO=ON \
-DTrilinos_ENABLE_Sacado=ON  \
-DTrilinos_ENABLE_Teuchos=ON \
-DTrilinos_ENABLE_MueLu=ON   \
-DTrilinos_ENABLE_ML=ON  \
-DTrilinos_VERBOSE_CONFIGURE=OFF \
-DTPL_ENABLE_MPI=ON  \
-DBUILD_SHARED_LIBS=ON   \
-DCMAKE_VERBOSE_MAKEFILE=OFF \
-DCMAKE_BUILD_TYPE=RELEASE   \
-DCMAKE_INSTALL_PREFIX:PATH=$HOME/share/trilinos \
../

make install


I got error message as follows...




"

The Fortran compiler identification is unknown

CMake Error at 
cmake/tribits/core/package_arch/TribitsGlobalMacros.cmake:1638 
(ENABLE_LANGUAGE):

  No CMAKE_Fortran_COMPILER could be found.


  Tell CMake where to find the compiler by setting either the environment

  variable "FC" or the CMake cache entry CMAKE_Fortran_COMPILER to the full

  path to the compiler, or to the compiler name if it is in the PATH.

Call Stack (most recent call first):

  cmake/tribits/core/package_arch/TribitsProjectImpl.cmake:188 
(TRIBITS_SETUP_ENV)

  cmake/tribits/core/package_arch/TribitsProject.cmake:93 
(TRIBITS_PROJECT_IMPL)

  CMakeLists.txt:89 (TRIBITS_PROJECT)

"


What's the problem? 

and I would very thank if someone tell me how to change deal.ii 
configuration to set on Trilinos after that


Thank you! 


-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Some questions on Nonlinear Problems

[JAEKWANG KIM]

Hi, all. 
I got a question on the symmetry of system matrices. 

I am looking for ways to solve Non-linear PDE using deal.ii and Dr. 
Bangerth's video lecture from 31.5~31.7 is greatly helpful to understand 
basic of this part. 

In many elementary codes in tutorial lists, we have used Conjugate Gradient 
Method to solve the matrix. 
CG Method basically assumes that matrix A is symmetric. 
I also checked that deal.ii offers other solver method, GMRES solver for 
non-symmetric system matrix A. 

Question 1. My question is where will the non-symmetry arise? 

To be more specific, let's take a example of "linearized" non-linear 
minimal surface equation with Picard iteration (as Dr. Bangerth's lecture 
31.65)

<https://lh3.googleusercontent.com/-TYzXkQCHcTA/WAksZOt6tjI/A7E/kX0e1wXWcfMcKc-dDYeq2rOOmXi-YtHjwCLcB/s1600/Screen%2BShot%2B2016-10-20%2Bat%2B3.42.54%2BPM.png>
where u is scalar function.  
Because we use previous solution u_k , the term "1/sqrt(1+(nabla_u)^2)" 
will work as a scalar coefficient that varies over the domain. 

Does this fact will make our system matrix non-symmetric? 

I have tried some mathematics to this.and I think that this might be 
non-symmetric , but still cannot sure on this...
(and think that if it is non-symmetric, than we would better to use GMRES 
Solver instead of CG Solver)


Question 2. Can we expect our numerical error will also behave ~O(h^p) for 
non-linear system that includes iterative method?

Not only mesh spacing and shape function approximation degree, but also 
many other parameters might influence error. 
(e.g. mapping degree, and degree of gauss quadrature...)
For Non-linear problem with iterative method, the tolerance between 
previous solution and solution might also affect error  
However, if I use enough small tolerance between abs(u_{k}-u_{k+1}), will 
the error will also show ~O(h^p)?


Always thank you all and I am learning a lot form this amazing community!!

Jaekwang Kim  


   


-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Easy way to calculate second-deviatoric tensor

[JAEKWANG KIM]

Thank you all guys who replied here!! 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Easy way to calculate second-deviatoric tensor

[JAEKWANG KIM]

Hi, all. 

I am really receiving a lot of help from this group whenever I visit here. 

I have a question today. 


I am developing FEM code for analyzing non-newtonian fluids using deal.ii 
library, which is amazing tool. 

In many cases, I have to calculate 'The second invariant of rate-of-strain 
tensor'. 

I wonder there is a easy way to get this value. 

As in step-20 or other fluid problems, we get rate of strain tensor from 
the rank-2 tensor of symmetrize gradient. 


For example


std::vector<SymmetricTensor<2,dim> > symgrad_phi_u (dofs_per_cell);

.

symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient (k, q)


>From this symmetric gradient tensor, do we have one easy way to calculate 
the following "the second invariant of rate-of-strain tensor"? 

<https://lh3.googleusercontent.com/-k8SKQzE3aEk/WAQQnorR5lI/A60/4XJBD25ZRSAa2IYsdK4NCMOuY3QXJydmwCLcB/s1600/Screen%2BShot%2B2016-10-16%2Bat%2B6.42.40%2BPM.png>



Regards, 

Jaekwang Kim 
 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Relation between Solution Error Behavior and Polynomial Approximation Degree

[JAEKWANG KIM]

<https://lh3.googleusercontent.com/-dkfwJvgPEmc/V-00Y-c4tzI/A5s/pb0Z-IBRe-An_S39dSDzxWg1k5GJlEbKQCLcB/s1600/step7plotting.jpg>

Hi all, I have question on error behavior of FEM. 

I thought that the order of error is O(h^p) where h is a mesh-size and p is 
polynomial degree we use in approximation. 

So, I thought that if I plot an error with number of mesh in log-log scale, 
than the graph will show -p slope. 
However, I the error behaves little bit different from my expectation.

For example, I use a step7 tutorial program (which solves Helmholtz 
decomposition and compares the FEM solution with exact solution.) 

The error curve showed more steep slope whenever I increase polynomial 
degree approximation however, the slope is not (-p). 
I reached slope (-3) when I used fifth-degree polynomial approximation...   
You can check this behavior in attached picture. 

Until now, I have considered, 

1. Mapping(From reference cell to real cell) degree (which is originally 
set to 1 but I used higher mapping) 
2. Instead of Qgauss quadrature, I am using QgaussLobatto Quadrature for 
any integration over cells. 
3. Shape function , again I tried to use QgaussLobatto node point for 
this) 

is there any suggestion that I need to fix more? 
or my first prediction that the slope will show '-p' or error will just 
behave O(h^p) was wrong?

I am always thank you for all guys!

Jaekwang Kim  

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Gauss-Lobatto node distribution on FESYSTEM?

[JAEKWANG KIM]

I want to use Gauss-Lobatto node on my fesystem for vector valued problems. 

While seeing step-48, I see that how Gauss-Lobatto node is used but I 
couldn't figure out how I should use it for vector value problems. 

To be specific, 

In step 48, fe is declared in main class as follow 


  template 

  class SineGordonProblem

  {

  public:

SineGordonProblem ();

void run ();


  private:

ConditionalOStream pcout;


void make_grid_and_dofs ();

void output_results (const unsigned int timestep_number);


#ifdef DEAL_II_WITH_P4EST

parallel::distributed::Triangulation   triangulation;

#else

Triangulation   triangulation;

#endif

*FE_Qfe;*

DoFHandler  dof_handler;

ConstraintMatrix constraints;

IndexSet locally_relevant_dofs;

 

  };

and than later, the constructor specified it as Gauss Lobbato point. 


  template 

  SineGordonProblem::SineGordonProblem ()

:

pcout (std::cout,

   Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0),

#ifdef DEAL_II_WITH_P4EST

triangulation (MPI_COMM_WORLD),

#endif

*fe (QGaussLobatto<**1>(fe_degree+1)),*

dof_handler (triangulation),

n_global_refinements (10-2*dim),

time (-10),

final_time (10),

cfl_number (.1/fe_degree),

output_timestep_skip (200)
  {} 


However, in my case, because I am using FESYSTEM 

it is declared in main class as follow... 

*FESystemfe; *


and constructor mentions taylor hood space as follow 


  StokesProblem::StokesProblem (const unsigned int degree)

:

degree (degree),

   * fe (FE_Q(degree+**1), dim,*

*FE_Q(degree), 1),*

dof_handler (triangulation)

  {}




I wonder how specific lines of coding makes my fesystem to use GaussLobatto 
node ... 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] VectorTools::integrate_difference function

[JAEKWANG KIM]

Thank you for all the repliers!
All of you gave me a good advice that I may consider deeply when I write my 
code from now
I have figured out what was wrong. 
Seeing at Step 20 code, I realized that I have to select component of 
solution before going to compare the difference. 

const ComponentSelectFunction

pressure_mask (dim, dim+1);

const ComponentSelectFunction

velocity_mask(std::make_pair(0, dim), dim+1);



Solution exact_solution;

Vector cellwise_errors (triangulation.n_active_cells());

   

QTrapez<1> q_trapez;

QIterated quadrature (q_trapez, degree+2);



// With this, we can then let the library compute the errors and 
output

// them to the screen:

VectorTools::integrate_difference (dof_handler, solution, 
exact_solution,

   cellwise_errors, quadrature,

   VectorTools::L2_norm,

   _mask);

const double p_l2_error = cellwise_errors.l2_norm();



VectorTools::integrate_difference (dof_handler, solution, 
exact_solution,

   cellwise_errors, quadrature,

   VectorTools::L2_norm,

   _mask);

const double u_l2_error = cellwise_errors.l2_norm();




-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Can we say that the higher order method, the more accurate?

[JAEKWANG KIM]

2016년 9월 18일 일요일 오전 10시 57분 16초 UTC-5, JAEKWANG KIM 님의 말:
>
>
> Hello, I am a starter of dealii and am learning a lot these days with the 
> help of video lectures and tutorial examples. 
>
> I modified step-22 code (stokes flow code) into my own problem, the flow 
> around sphere.
>
> and I intend to evaluate the drag force (which is analytically given by 
> stokes equation) 
>
> My code reached quite close to the value since the absolute error  : 
> abs(drag_calculated-drag_exact)/drag_exact is around 10^(-3)
>
> However, I expected that if I input higher 'degree' I will receive more 
> accurate result, but it didn't
>
> Obviously Q2 is better than Q1. and Q3 is better than Q2. But Q4 or Q4 is 
> not better than Q2 or Q3? 
>
> Is there any reason on this? 
>
> (To be specific, if i say degree 2 , that mean I use (2+1) for velocity, 
> (2) for pressure, and (2+2) for Gauss integral
>
>
> Thank you 
>
> Jaekwang Kim  
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Can we say that the higher order method, the more accurate?

[JAEKWANG KIM]

Thank you for all repliers!!

I tried fix code as you mentioned please see redline...

 template 

  void StokesProblem::assemble_system ()

  {

system_matrix=0;

system_rhs=0;


*QGaussLobatto<**2>   quadrature_formula(degree+2);*


FEValues fe_values (fe, quadrature_formula,

 update_values|

 update_quadrature_points  |

 update_JxW_values |

 update_gradients);



but I received an error message that







An error occurred in line <1831> of file 
 in 
function

void dealii::BlockMatrixBase::add(const 
size_type, const size_type, const value_type) [MatrixType = 
dealii::SparseMatrix]

The violated condition was: 

dealii::numbers::is_finite(value)

The name and call sequence of the exception was:

ExcNumberNotFinite(std::complex(value))

Additional Information: 

In a significant number of places, deal.II checks that some intermediate 
value is a finite number (as opposed to plus or minus infinity, or NaN/Not 
a Number). In the current function, we encountered a number that is not 
finite (its value is (nan,0) and therefore violates the current assertion.


This may be due to the fact that some operation in this function created 
such a value, or because one of the arguments you passed to the function 
already had this value from some previous operation. In the latter case, 
this function only triggered the error but may not actually be responsible 
for the computation of the number that is not finite.


There are two common cases where this situation happens. First, your code 
(or something in deal.II) divides by zero in a place where this should not 
happen. Or, you are trying to solve a linear system with an unsuitable 
solver (such as an indefinite or non-symmetric linear system using a 
Conjugate Gradient solver); such attempts oftentimes yield an operation 
somewhere that tries to divide by zero or take the square root of a 
negative value.


In any case, when trying to find the source of the error, recall that the 
location where you are getting this error is simply the first place in the 
program where there is a check that a number (e.g., an element of a 
solution vector) is in fact finite, but that the actual error that computed 
the number may have happened far earlier. To find this location, you may 
want to add checks for finiteness in places of your program visited before 
the place where this error is produced.One way to check for finiteness is 
to use the 'AssertIsFinite' macro.

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] VectorTools::integrate_difference function

[JAEKWANG KIM]

yes...I tried that before, but I get an error message as follow.. 

As you mentioned, I first learn this function in step 7 tutorial. 
But the problem is that I want to integrate my vector-valued-solution to 
exact-vector-valued solution, while step 7 compares scalar-solution. 

Or is there any ways to compare each solution component separately? as it 
did in step-7. 

I am really thank you for your advice!


An error occurred in line <235> of file 
 in 
function

static types::global_dof_index 
dealii::internal::DoFHandler::Policy::Implementation::distribute_dofs(const 
types::global_dof_index, const types::subdomain_id, DoFHandler &) [dim = 2, spacedim = 2]

The violated condition was: 

tria.n_levels() > 0

The name and call sequence of the exception was:

ExcMessage("Empty triangulation")

Additional Information: 

Empty triangulation


Stacktrace:

---

#0  2   libdeal_II.g.8.4.1.dylib0x000104ef307c 
_ZN6dealii8internal10DoFHandler6Policy14Implementation15distribute_dofsILi2ELi2EEEjjjRNS_10DoFHandlerIXT_EXT0_EEE
 
+ 300: 2   libdeal_II.g.8.4.1.dylib0x000104ef307c 
_ZN6dealii8internal10DoFHandler6Policy14Implementation15distribute_dofsILi2ELi2EEEjjjRNS_10DoFHandlerIXT_EXT0_EEE
 

#1  3   libdeal_II.g.8.4.1.dylib0x000104ef2d97 
_ZNK6dealii8internal10DoFHandler6Policy10SequentialILi2ELi2EE15distribute_dofsERNS_10DoFHandlerILi2ELi2EEERNS1_11NumberCacheE
 
+ 39: 3   libdeal_II.g.8.4.1.dylib0x000104ef2d97 
_ZNK6dealii8internal10DoFHandler6Policy10SequentialILi2ELi2EE15distribute_dofsERNS_10DoFHandlerILi2ELi2EEERNS1_11NumberCacheE
 

#2  4   libdeal_II.g.8.4.1.dylib0x000104eb6797 
_ZN6dealii10DoFHandlerILi2ELi2EE15distribute_dofsERKNS_13FiniteElementILi2ELi2EEE
 
+ 135: 4   libdeal_II.g.8.4.1.dylib0x000104eb6797 
_ZN6dealii10DoFHandlerILi2ELi2EE15distribute_dofsERKNS_13FiniteElementILi2ELi2EEE
 

#3  5   mystokes0x000101350e7d 
_ZN8MyStokes13StokesProblemILi2EE10setup_dofsEv + 397: 5   mystokes
0x000101350e7d 
_ZN8MyStokes13StokesProblemILi2EE10setup_dofsEv 

#4  6   mystokes0x000101341a52 
_ZN8MyStokes13StokesProblemILi2EE3runEv + 386: 6   mystokes
0x000101341a52 _ZN8MyStokes13StokesProblemILi2EE3runEv 

#5  7   mystokes0x000101340aa5 main + 69: 7 
  mystokes0x000101340aa5 main 

#6  8   libdyld.dylib   0x7fff95b575ad start + 1: 8 
  libdyld.dylib   0x7fff95b575ad start 

#7  9   ??? 0x0001 0x0 + 1: 9   
??? 0x0001 0x0 

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] VectorTools::integrate_difference function

[JAEKWANG KIM]

thanks!! Full error message is as follow 


*/Users/jaekwangjk/repos/trial/sphere-deal.ii/mystokes.cc:1258:7: **error: *
*no*

*  matching function for call to 'integrate_difference'*

  VectorTools::integrate_difference (dof_handler,

*  ^*

*/Users/jaekwangjk/repos/trial/sphere-deal.ii/mystokes.cc:1220:7: note: *in

  instantiation of member function

  'MyStokes::StokesProblem<2>::process_solution' requested here

  process_solution ();

*  ^*

*/Users/jaekwangjk/repos/trial/sphere-deal.ii/mystokes.cc:1298:21: note: *in

  instantiation of member function 'MyStokes::StokesProblem<2>::run'

  requested here

   flow_problem.run ();

*^*

*/Users/jaekwangjk/Programs/dealii/include/deal.II/numerics/vector_tools.h:1940:8:
 
note: *

  candidate template ignored: could not match 'Function' against 
'Solution'

  void integrate_difference (const DoFHandler<dim,spacedim> ,

*   ^*

*/Users/jaekwangjk/Programs/dealii/include/deal.II/numerics/vector_tools.h:1968:8:
 
note: *

  candidate template ignored: could not match 'dealii::hp::DoFHandler'

  against 'dealii::DoFHandler'

  void integrate_difference (const hp::DoFHandler<dim,spacedim> ,

*   ^*

*/Users/jaekwangjk/Programs/dealii/include/deal.II/numerics/vector_tools.h:1925:8:
 
note: *

  candidate function template not viable: requires at least 7 
arguments, but

  6 were provided

  void integrate_difference (const Mapping<dim,spacedim>,

*   ^*

*/Users/jaekwangjk/Programs/dealii/include/deal.II/numerics/vector_tools.h:1953:8:
 
note: *

  candidate function template not viable: requires at least 7 
arguments, but

  6 were provided

  void integrate_difference (const hp::MappingCollection<dim,spacedim> 
,




I will check the mapping as you advised for my other questions.!!

2016년 9월 20일 화요일 오전 9시 51분 13초 UTC-5, Martin Kronbichler 님의 말:
>
> Jaekwang,
>
> can you post the full error message? To me the call looks reasonable (this 
> is a copy from step-7) and the function definition looks good as well.
>
> Ah, and in case you want to use this in the context of your other question 
> regarding high order methods: Put a "mapping" as first argument to this 
> call and make sure to use QGauss(fe.degree+2) or so, which increases 
> the order of quadrature as the element order gets higher.
>
> Best,
> Martin
>
> On 09/20/2016 04:46 PM, JAEKWANG KIM wrote:
>
>  
>
> I wanted to evaluate Vectortool::integrate_difference to estimate my 
> error. 
>
>
>
> but I couldn't figure out what I need to enter for the third component. (I 
> marked it as red) 
>
>
>
>
>  VectorTools::integrate_difference (dof_handler,
>
>  solution,
>
>  *Solution<dim+**1**>()*,
>
>  difference_per_cell,
>
>  QGauss(3),
>
>  VectorTools::L2_norm);
>
>
> I have solution vector over my domain,(by the way I am solving vector 
> valued problems...)  
>
> I want to compare this with continuous exact function, "May be 
> Solution<dim+1>" , I have generated it as follow. 
>
> However, I got a error message that 
>
>
> */Users/jaekwangjk/repos/trial/sphere-deal.ii/mystokes.cc:1258:7: **error: 
> **no*
>
> *  matching function for call to 'integrate_difference'*
>
>
> What's might be wrong with my code..
>
>
>  template 
>
> class Solution : public Function
>
> {
>
> public:
>
> 
>
> Solution () : Function(dim+1) {}
>
> 
>
> //const unsigned int b_id;
>
> virtual double value (const Point   ,
>
>   const unsigned int component = 0) const;
>
> 
>
> virtual void vector_value (const Point ,
>
>Vector   ) const;
>
> 
>
> };
>
> 
>
> 
>
> template 
>
> double
>
> Solution::value (const Point  ,
>
> const unsigned int component) const
>
> {
>
> Assert (component < this->n_components,
>
> ExcIndexRange (component, 0, this->n_components));
>
> 
>
> double U=-1.; //inflow velocity
>
> double obj_rad=0.2; //Object Radius
>
> double r_sph = sqrt( pow(p[0],2) + pow(p[1],2) );
>
> double phi = (-1) * acos(p[1]/r_sph);
>
> 
>
> double q_r_sph=U*cos(phi)*( 1

Re: [deal.II] Re: Can we say that the higher order method, the more accurate?

[JAEKWANG KIM]

Can I ask more about "the limits of our implementation of our quadrature 
formulas?". 
I wonder when it usually happens. 


Once I calculated drag coefficient, with Q1, my error is 10% compare to 
exact solution. 
However, I can 1.7% error when I use Q2 which is significant decrease!
At Q3, I get 0.4% and At Q4 it starts to increase again 0.5%. 
If I go higher, then my error is more than 100%. I really want to 
figure out why this happens

To summarize 
>From Q1~Q3... it shows significant decrease in error
but it is not anymore at Q4 and Q5 

2016년 9월 19일 월요일 오전 8시 8분 21초 UTC-5, Guido Kanschat 님의 말:
>
> You will also encounter the limits of our implementation of our quadrature 
> formulas and our polynomial evaluation. Both are not stable for high order, 
> and thus at some point approximation will stall.
>
> Best,
> Guido
>
>
>
> -- 
> Prof. Dr. Guido Kanschat
> Interdisziplinäres Zentrum für Wissenschaftliches Rechnen
> Universität Heidelberg
> Im Neuenheimer Feld 368, 69120 Heidelberg
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Can we say that the higher order method, the more accurate?

[JAEKWANG KIM]

I also wonder whether their is a limitation in higher order method that 
deal offer. 

If no, then we can use Q10 method, just putting degree value as 10?  

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Can we say that the higher order method, the more accurate?

[JAEKWANG KIM]

This is how my mesh looks like. 
I am using half domain of 2D sphere assuming that it would symmetrical 
along phi domain. 
Tho my domain is 2d, I am describing 3D phenomena. I am using cylindrical 
coordinate. 
So... in real horizontal line mean r-axis and vertical line is z-axis. 

When r=0 (so the left straight line), I used boundary condition that u_r=0 
(because there will be no flow across this axis if flow is symmetric) 

other boundary, sphere boundary condition is 0 velocity, and top,end,right 
side boundary condition is given by analytical solution (it was expressed 
in spherical coordinate, so I transferred it into cylindrical coordinate) 

To, Daniel Arndt. 
Also, I didn't understand what you meant by 
*"Are you using a Mapping of the same order as your velocity ansatz space 
is?"*


Thanks for reply!!!

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.