[deal.II] Re: Non Homogenous Boundary Conditions on a Cylinder

2020-06-30 Thread Simon Sticko
Hi,
>From the figure you can't really tell whether the boundary condition is 
fulfilled or not. If you want to check if the normal derivative is 
reasonable (almost equal to 0) at the boundary, the easiest is probably to 
plot the solution on a line in the radial direction. This can be done by 
using
Filters -> Data Analysis -> Plot Over Line in Paraview.

I don't quite understand what you mean when you say that you have tried 
changing the value of component. Component is intended to be used for 
vector valued problems. If you are solving Poisson, it doesn't make sense 
to set it to something different than 0. 

Note that, right now you are calling the function 
NeumannBoundaryCondition::value with a single argument (of type Point<3>). 
This means that you get the default argument for component: component=0.

Best,
Simon

On Tuesday, June 30, 2020 at 2:59:20 AM UTC+2, Samuel Rodriguez wrote:
>
> Hi Simon,
>
> It seems I might not have understood my own question. You answered the 
> question I was trying to ask. Would you be able to help me out with one 
> more thing? I followed your advice but the boundary conditions do not seem 
> bedo not seem to be getting the correct boundary conditions, so I don't 
> think they are being applied correctly. Here is the simple piece of code 
> that I wrote:
>
> class NeumannBoundaryCondition : public Function<3>
> {
> public:
>   double
>   value(const Point<3> ,  const unsigned int component = 0) const
>   {
> return component;
>   }
> };
>
>
> This neumann boundary condition is only meant to apply a constant value of 
> zero on the sides. Here is how I implemented it:
>
> QGauss<3> quadrature_formula(fe.degree + 1);
> QGauss<3 - 1> face_quadrature_formula(fe.degree + 1);
> FEFaceValues<3> fe_face_values(fe,
>face_quadrature_formula,
>update_values | update_quadrature_points |
>update_normal_vectors |
>update_JxW_values);const NeumannBoundaryCondition 
> boundary_condition;
>
> const unsigned int n_face_q_points = face_quadrature_formula.size();
>
>
> for(unsigned int face_number = 0; face_number < GeometryInfo<3>::
> faces_per_cell; ++face_number)
> {
>   if(cell->face(face_number)->at_boundary() && (cell->face(face_number)->
> boundary_id() == 0))
>   {
>fe_face_values.reinit(cell, face_number);
>for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) {
>  const double neumannvalue = boundary_condition.value(fe_face_values.
> quadrature_point(q_point));
>  for (unsigned int i = 0; i < dofs_per_cell; ++i) {
>const unsigned int component_i = fe.system_to_component_index(i).
> first;
>cell_rhs(i) += (neumannvalue *
>fe_face_values.shape_value(i, q_point) *
>fe_face_values.JxW(q_point));
>}
>  }
>  
>   }
> }
>
>
> This seems correct to me, but the code is showing a cyldiner that looks 
> like this:
>
>
> [image: Screen Shot 2020-06-29 at 5.50.19 PM.png]
> which doesn't appear to be the right solution. Even when I change the 
> value of component to something ridiculous such as 10,000,000 the solution 
> on the cylinder does not change. Would you be able to tell me what I am 
> doing wrong in the implementation?
>
> - Samuel
>
>
>
> On Friday, June 26, 2020 at 2:21:19 AM UTC-7, Simon Sticko wrote:
>>
>> Hi,
>> I'm not sure I understand your question correctly, but if you want to 
>> impose an inhomogeneous Neumann boundary condition you would typically 
>> write a small class representing the boundary condition you want:
>>
>> class NeumannBoundaryCondition : public Function<3>
>> {
>> public:
>>   double
>>   value(const Point<3> , const unsigned int component = 0) const
>>   {
>> // Implement the boundary condition you want here.
>> return 1;
>>   }
>> };
>>
>> Then just use this when you assemble:
>>
>> const NeumannBoundaryCondition boundary_condition;
>>
>> const Point<3>  = fe_face_values.quadrature_point(q_point);
>>
>> const double neumann_value = boundary_condition.value(point);
>>
>> cell_rhs(i) += (neumann_value *
>> fe_face_values.shape_value(i, q_point) *
>> fe_face_values.JxW(q_point));
>>
>>
>> Does that answer your question?
>>
>> Best,
>> Simon
>>
>> On Friday, June 26, 2020 at 6:02:48 AM UTC+2, Samuel Rodriguez wrote:
>>>
>>> Hello Everyone,
>>>
>>> I am a masters student barely getting into working with the very heavy 
>>> theoretical foundations of using deal.ii. I have been working on a research 
>>> project that uses a code that creates a cylinder with Dirichlet boundary 
>>> conditions on the hull/bottom faces and a Neumann boundary condition on the 
>>> top face Recently I have been trying to change the Dirichlet boundary 
>>> condition on the hull of the cylinder to a Neumann boundary condition. Our 
>>> code can be found here inside of QuasistaticBrownianThermalNoise.cpp:
>>>
>>> Numerical Coating Thermal Noise 
>>> 
>>>
>>> I 

[deal.II] Is there an advantage/disadvantage/difference to using FEInterfaceValues (tutorial step 12) instead of the MeshWorker framework (tutorial step 12b) when it comes to DG?

2020-06-30 Thread CopyCat
Hello,

I am in the process of learning DG and dealii. 
I would like to implement a DG solver that would solve the 
convection-diffusion equation and later implement another code to solve the 
incompressible NS equation using DG as well. Tutorials 12 and 12b use 
different approaches to solve the pure convection equation, and so I was 
wondering which of these two approaches is better suited for my 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/c127437d-778e-4a1a-a5b0-72b5c4c4b51co%40googlegroups.com.


Re: [deal.II] GMesh 1D mesh input error

2020-06-30 Thread Wolfgang Bangerth

On 6/26/20 11:05 AM, Victoria W. wrote:
This fixed my problem in 2d, but I'm still having issues in 3d.  Any other 
suggestions for getting a .msh file read in when it's a 2d mesh in 3d space?  
Posting the mesh I want to use here, with the type 15 elements removed.  I've 
also tried other software and file formats, running into the same issue 
previously posted about .unv files produced with Salome, so I appreciate any 
suggestions on reading in a mesh!


Victoria,
I assume that you are talking about the following error:

 

An error occurred in line <2674> of file 
 in function
static void 
dealii::internal::TriangulationImplementation::Implementation::create_triangulation(const 
std::vector >&, const std::vector >&, 
const dealii::SubCellData&, dealii::Triangulation<2, spacedim>&) [with int 
spacedim = 3]

The violated condition was:
line->boundary_id() != numbers::internal_face_boundary_id
Additional information:
The input data for creating a triangulation contained information about a 
line with indices 9 and 33 that is described to have boundary indicator 0. 
However, this is an internal line not located on the boundary. You cannot 
assign a boundary indicator to it.


If this happened at a place where you call 
Triangulation::create_triangulation() yourself, you need to check the 
SubCellData object you pass to this function.


If this happened in a place where you are reading a mesh from a file, then you 
need to investigate why such a line ended up in the input file. A typical case 
is a geometry that consisted of multiple parts and for which the mesh 
generator program assumes that the interface between two parts is a boundary 
when that isn't supposed to be the case, or where the mesh generator simply 
assigns 'geometry indicators' to lines at the perimeter of a part that are not 
supposed to be interpreted as 'boundary indicators'.




The issue here is similar to the previous one, just that now you have a *line* 
in the interior of a 2d mesh (where previously you had a vertex interior to a 
1d mesh) for which the input file specifies boundary values -- even though 
that line is not actually at the boundary. These lines are now "type 1" 
entities in the language of the GMSH and look like this:

  33 1 0 1 2 10 34
(GMSH numbers vertices starting at 1, whereas deal.II uses 0, so this is the 
line that the error message above complains about: the one with vertices 9 and 
33.)


In your file, there are 208 such lines (listed under the "$ELM" line, and 
numbered 33 to 240). If you remove these lines and adjust the number under 
$ELM from 802 to 594, the file can be read successfully and looks like the one 
attached.


Best
 W.

--

Wolfgang Bangerth  email: bange...@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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/78882fe6-ff7e-87e1-0c46-a40c25b50a86%40colostate.edu.


Re: [deal.II] Is there an advantage/disadvantage/difference to using FEInterfaceValues (tutorial step 12) instead of the MeshWorker framework (tutorial step 12b) when it comes to DG?

2020-06-30 Thread Wolfgang Bangerth



I would like to implement a DG solver that would solve the 
convection-diffusion equation and later implement another code to solve the 
incompressible NS equation using DG as well. Tutorials 12 and 12b use 
different approaches to solve the pure convection equation, and so I was 
wondering which of these two approaches is better suited for my problem?
The two approaches were developed for different purposes: FEInterfaceValues as 
an approach to compute jump terms for DG methods in the same way as use 
FEValues and FEFaceValues for evaluating the shape functions on cells and 
faces. On the other hand, MeshWorker was meant as a framework that would make 
the assembly of linear and bilinear forms over cells and interfaces easier. If 
FEInterfaceValues had been around at the time MeshWorker was written, it would 
have been one of the building blocks used there.


So there is no real equivalency between the two. The equivalence is more 
between MeshWorker framework and the MeshWorker::mesh_loop() method used in 
step-12 and step-47: As has been expressed many times on this mailing list, 
MeshWorker has a good underlying idea, but it is not well documented or tested 
and none of the people who frequently answer on the mailing list know it well. 
In other words, its higher level functionality is in essence unmaintained. On 
the other hand, MeshWorker::mesh_loop() *is* well documented and supported, 
and has an interface that is relatively widely understood here. My preference 
would have been to just remove step-12b when the previous version was 
rewritten into what step-12 is now, rather than renaming it step-12b. (And do 
the same for step-16b.) I consider step-12b and step-16b as obsolete.


Best
 W.

--

Wolfgang Bangerth  email: bange...@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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/cf5ffbb2-485e-7deb-f5f7-8633af2c1c35%40colostate.edu.